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I. INTRODUCTION 


This final report is for work carried out under Grant No. NCC2-5072 during 
the period from August 1, 1994 to September 30, 1996. 

The objective of this study was to modify an existing parallel particle code, 
based on the direct simulation Monte Carlo (DSMC) method, to include a Navier- 
Stokes (NS) calculation so that a hybrid solution could be developed. The former was 
developed by McDonald, Fallavollita and Dahlby at Stanford University for use on 
the Intel Gamma, Delta, and Paragon computers. The need for a hybrid DSMC/NS 
solution arises, for example, in situations in which a cold blunt body is placed in a 
high enthalpy flow. The very high stagnation temperature and pressure produced, 
together with the cold body, cause the gas density near the surface to become quite 
high; the density ratio across the thermal layer can be 30 or greater. Because of 
the high density and the large gradients present, use of the DSMC method alone 
becomes computationally costly. The aim was to make use of a Navier-Stokes code 
to handle all near-continuum sub-regions, such as the thermal layer or a boundary 
layer in general. 

In carrying out the work, it was determined that the following five issues had to 
be addressed before extensive program development for 3D capability was pursued: 

( i ) find a set of one-sided kinetic fluxes that are fully compatible with the DSMC 
method, so that they may be used in interfacing the two solution schemes; 

( ii ) develop a finite volume scheme that makes use of these one-sided kinetic fluxes; 

(in) make use of the one-sided kinetic fluxes together with DSMC type boundary 
conditions at a material surface so that velocity slip and temperature slip arise 
naturally for near-continuum conditions; 

(iv) find a suitable sampling scheme so that the values of the one-sided fluxes pre- 
dicted by the NS solution at an interface between the two domains can be 
converted into the correct distribution of particles to be introduced into the 
DSMC domain; and 

( v ) carry out a suitable number of tests to confirm that the developed concepts are 
valid, individually and in concert for a hybrid scheme. 
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The first three issues were addressed by S. Chou and elegantly handled in his 
Ph.D. thesis [1], and in a follow-on paper which was recently accepted for publication 
in the Journal for Computational Physics [2]. It became possible to address the fourth 
issue directly by making use of the analytic expressions developed by Chou, and its 
implementation is discussed in the two theses by D. Dahlby [3] and by T. Lou [4]. The 
fifth issue was studied by Lou [4] and by Lou et al. [5] for one-dimensional geometries, 
and by Dahlby [3] and Duttweiler [7] for two dimensional geometries. Presentations 
of the work were made by Dahlby [6] and most recently by Duttweiler [7]. All of the 
listed work was supported, to varying degrees, by the grant. 

Because much of the work has been reported in student theses and represents 
many pages of text, only a short review of the relevant topics, which can be used as 
a guide to the publications themselves, will be given. On the other hand, the two 
papers [2] and [5] have yet to appear as archive journal publications and so copies are 
attached as Appendices A and B. 


II. Split Kinetic Fluxes, Numerical Studies 

Starting with the kinetic-theory definition for the split kinetic fluxes, their 
derivation and development for the Chapman-Enskog velocity distribution function 
were carried out by Chou [1], and a listing of the resultant set of relations for the 
split fluxes is given as Eqs. (38)-(44) in Appendix A. The appropriate boundary condi- 
tions for the split fluxes are described by Eqs. (62)-(68) in the same appendix. These 
boundary conditions are flux-type boundary conditions and therefore lead to both 
temperature slip and velocity slip at a material surface. A version of the monotone 
upstream-centered scheme for conservation laws (MUSCL) was used in applying the 
split kinetic flux relations to obtain numerical solutions for the NS equations for var- 
ious test flows, as described on page 11 of Appendix A. Numerical tests on the split 
kinetic fluxes were carried out for the highly viscous and heat conducting flow found 
in a normal shock wave profile (App. A, Figs. 5 and 6) and for the opposite limit 
represented by the Euler flow in the shock tube problem (Figs. 7 and 8). As shown 
in the figures, comparisons were made with numerical solutions to the Navier-Stokes 
equations based on the established methods of Roe and Jameson (SLIP2). Addition- 
ally, tests involving a solid surface were used to study the thermal boundary layer 
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found near an impulsively started piston (App. A, Figs. 1 and 9), which leads to 
temperature slip, along with tests involving the viscous, heat conducting flow found 
in the slider-plate or lid-driven cavity problem (Figs. 10 - 12), which introduces both 
temperature and velocity slip. As seen in the figures, all comparisons are uniformly 
good. 


III. Confirmation of Split Kinetic Fluxes and Slip 

Because the split kinetic fluxes are derived for the nonequilibrium state repre- 
sented by the Chapman- Enskog velocity distribution function, it is expected that the 
values predicted by the analytic expressions should agree with results of simulations 
using the direct simulation Monte Carlo (DSMC) method, as long as the Knudsen 
number, velocity gradients, and temperature gradients are suitably small in the sim- 
ulated flow. In addition, because the boundary conditions used in the two methods 
are physically identical, one would expect good agreement for all slip values as well. 
Studies of this hypothesis were carried out by Dahlby [3] and by Lou [4], and a paper 
summarizing the work [5] is reproduced as Appendix B. The tests were conducted for 
conditions of flow (Knudsen number and dimensionless time) in which both methods 
(DSMC and NS) alone were capable of providing the physically correct prediction, 
and thus should agree. Three test cases were studied. 

The first case considered an impulsively started piston in a stationary equilib- 
rium gas, and confirmation that the temperature and density profiles in the thermal 
layer agree after roughly 10 collision times is shown in Figs. 2 and 3 in Appendix 
B. Because the numerical value of temperature slip for this case is small and on the 
order of 10%, and because sampling in the DSMC method is severely limited in a 
nonsteady problem where time averaging in not permitted, accurate DSMC values 
are difficult to obtain. Still, the comparison shown in Fig. 4 of Appendix B is re- 
markably good. The split fluxes for mass, momentum and energy were recorded at 
the wall in the DSMC method by monitoring the particles that cross from the gas 
into the wall region during a tinle step. Comparisons between theory and simulation 
are presented in Figs. 5 and 6, showing outstanding agreement for the split fluxes. 

The second test considered an impulsively started flat plate moving parallel to 
its surface in a stationary equilibrium gas. Here the nonequilibrium is caused by the 
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velocity field and to a lesser degree by the temperature field. The velocity slip at 
the plate surface is displayed in Fig. 8 of Appendix B for the two methods showing 
remarkably good agreement. Likewise, the split fluxes are given in Figs. 9 and 10 
offering equally good agreement. 

The third test considered the steady state solution for the slider plate problem, 
where a two-dimensional square cavity filled with a gas is closed by a moving plate. 
The most impressive comparisons are given by the tangential velocity, Fig. 13, and the 
split fluxes, Fig. 14, at the surface of the cavity. Together, these three tests verify that 
the analytic expressions we are using for the split kinetic fluxes are in full agreement 
with results of simulations using the direct simulation Monte Carlo method. 


IV. Sampling from the Chapman-Enskog Distribution 

To implement a DSMC/NS hybrid scheme based on the passing of split-fluxes, 
one must be able to use the data from the NS solution (at the interface and at each 
time step) to create a proper sample of particles for insertion into the DSMC region. 
This collection of particles must carry the correct mass, momentum, and energy and 
in order to do so must represent a sample from a weighted Chapman-Enskog velocity 
distribution. A straightforward selection/rejection scheme is far too computationally 
costly, for a complex 3D function such as the CE distribution, and a less expensive 
method must be found. Conceptually, a superior procedure can be understood by a 
review of Fig. 1, where the equivalent 2D problem is considered. Taking Cx and Cy 
to be the two components of molecular velocity and assuming we are interested in the 
positive directed flux along the x-axis, then one need only consider the portion of the 
velocity distribution function shown in the first row. However, if we are interested 
in the distribution that makes up the one-way flux, then we need the product of Cx 
and f(Cx,Cy) which is the weighted distribution shown in the second row. If this 
distribution is integrated over the entire Cx domain, one then obtains a probability 
density that is a function of Cy alone, namely the marginal distribution, as shown 
in the dashed plane of the mesh plot of the third row. On creating the indefinite 
integral of this function, as displayed in the right-hand column, one obtains the 
probability that the y-component of the molecular velocity is smaller than some Cy 
for all values of Cx. This function can then be used to randomly pick a value of Cy 


4 



for a sample particle by first selecting a random number 0 < Rl < 1 and then finding 
the corresponding value of Cy from the plot. This value for Cy is then used in the 
distribution shown in the second row to obtain a function of C x alone, as depicted 
by the dashed plane in the mesh plot of the fourth row. 



Fig. 1 . Procedure for selecting a pair of random velocities Cx 
and Cy, from the weighted Chapman-Enskog distribution. 

On creating the indefinite integral of this function, as displayed in the right-hand 
column, one obtains the probability that the x-component of the molecular velocity 
is smaller than some Cx for the particular value of Cy selected. This function can 
then be used to randomly pick a value of Cx for the sample particle, by first selecting 
a random number 0 < R2 < 1 and then finding the corresponding value of Cx from 
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the plot. In summary, we pick two random numbers Rl and R2 to select two random 
velocities Cx and Cy for our first sample particle; and these steps are then repeated 
for all particles to be sampled. 

A typical sample for 200 particles is shown in Fig. 2, superposed on a contour 
plot of the weighted distribution itself, verifying the effectiveness of the procedure. 
We are able to use this procedure because all of the analytic expressions for the 
integrals needed, both definite and indefinite, were previously obtained by Chou in 
the process of developing his analytic results for the split kinetic fluxes. 



Fig. 2. A representative example for creating a sample of 200 
particles using the procedure outlined in Fig. 1. 


Although the method itself represents a general mathematical concept that is not new, 
the ability to carry it out for the Chapman-Enskog distribution is new because of the 
availability of Chou’s analytic results. Even though the sampling procedure outlined 
above is considerably more efficient than the simple selection/rejection scheme, it is 
still fairly intensive and steps are being taken to reduce the computational effort still 
further. Some of these steps were taken by Duttweiler in creating a very efficient C 
program to implement the algorithm. Duttweiler’s coding was used by Dahlby and 
by Lou in their respective thesis work. 
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V. Testing of a Flux-Based Hybrid Scheme 

Considerable testing of a flux-based hybrid scheme was carried out by Lou for 
two nonsteady one-dimensional problems and by Dahlby for a steady two-dimensional 
problem. The first problem studied by Lou was the impulsively started piston, where 
the interface for the hybrid was placed at a fixed position 12.5 mean free path lengths 
ahead of the piston surface. Some of his results (taken from his thesis [4]) for the case 
of a monatomic gas and M pis ton = 1-0 are shown in Fig. 3. 


PISTON, Time = 4.635 



Fig. 3. A hybrid solution for an impulsively started piston, 
with the piston surface on the right, and the interface (dotted 
vertical line) at 12.5 mean free path lengths ahead of the piston. 
Dimensionless time is based on the upstream collision time. 

The objective was to select conditions and a time scale for which one could observe 
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the wave front as it passed through the interface, to see if any imperfections appear 
at the point where the two solutions are joined. Far upstream and far downstream 
from the shock wave the conditions are essentially at equilibrium and no difficulty 
should be expected, but when the shock wave passes the interface location then severe 
nonequilibrium conditions are present, giving the matching algorithm its most severe 
test (intermediate time shown). Also shown in the figure, as a dashed curve, is a 
pure NS solution which is hardly distinguishable from the hybrid, indicating that 
the matching algorithm seems to be working rather well. For the later time shown, 
the hybrid arrangement is optimum for this problem, in that the shock wave profile 
itself can best be modeled by the DSMC method while the thermal layer at the 
piston surface can best be handled by the NS solution. The density profile for the 
intermediate time is quite similar in many respects to the density profile along the 
stagnation streamline in a supersonic flow past a blunt cold body, a situation of 
interest to us because it becomes a very useful application of the hybrid scheme. 

The second problem studied by Lou was the impulsively started flat plate, where 
the interface for the hybrid was placed at a fixed position 10 mean free path lengths 
above the plate surface. Some of his results (taken from his thesis [4]) lor the case of a 
monatomic gas and M p i a t e = 1.0 are shown in Fig. 4. Here again, the objective was to 
select conditions and a time scale for which one could observe the boundary layer as 
it grew and passed thought the interface, to see if any imperfections appeared at the 
point where the two solutions join. In this case the nonequilibrium was principally the 
result of stress with heat flux playing a lesser role. As can be seen, the matching again 
appears to be working very well. In this case it is also clear that the flux boundary 
conditions used in the KFVS method for the NS equations is accounting for velocity 
slip at the plate surface. The slip value is initially large and then diminishes as time 
progresses, reaching a level of about 10% for the latest time given. 

For both the impulsively started piston and the impulsively started flat plate, 
the conditions were chosen so that a pure NS solution and a pure DSMC solution 
would predict the same results. This is possible as long as the Knudsen number is 
sufficiently small and the nonequilibrium is not too large. Because of this, we were 
free to place the interface at most any position in the hybrid solution (positioning was 
not critical); and we were always in a position to use either pure solution as a reference 
in judging the hybrid. In Lou’s thesis [4] many of these comparisons are made, even 
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with a reversed DSMC/NS orientation. Also, in both cases it was necessary to study 
extremely short times in order for the nonequilibrium to be large, so that the hybrid 
scheme could be fully tested. Typically, a time of about 10 to 20 collision times was 
used, clearly a very short time in terms of physical processes. 


FLAT PLATE, Tune = 6.18 



U/Uplate 

Fig. 4. A hybrid solution for an impulsively started flat plate, 
with the plate surface at the bottom, and the interface (indi- 
cated by the dotted horizontal line) at 10 mean free path lengths 
above the plate. Dimensionless time is based on the freestream 
collision time. 

The third problem studied was the lid-driven cavity problem, or the slider plate 
problem. Large nonequilibrium develops in the two corners between the lid and the 
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cavity wall, and consequently, it is reasonable to assign the DSMC method a small 
strip alongside the lid, including the two corners; and then allow NS to handle the 
rest, where the degree of nonequilibrium is considerably smaller. The reason for 
this choice is that only DSMC alone can handle conditions where velocity slip and 
temperature slip may be large, and where the physical scale in the two corners may 
effectively lead to rarefied conditions. The hybrid solution for the pressure field is 
shown in Fig. 5, where the hybrid interface was placed at 50% of the depth of the 
cavity which can be visually located by the change in ’texture’ of the surface. 



Fig. 5. A hybrid solution for the lid-driven cavity problem 
showing the pressure field for Kn = 0.01, Mud = 1.0 and a 
monatomic gas. The shearing lid is on the far wall as indicated. 


Because Dahlby carried out this study [3] at a fairly early point in our work, before we 
fully developed the means for sampling from the Chapman-Enskog velocity distribu- 
tion function, he chose to place the interface at a position in the flow for which it was 
assured that nonequilibrium would be small, namely, the 50% location as shown. In 
this location, one is able to use the much simpler Maxwellian distribution, for which 
the proper sampling scheme is well known and straightforward. This test confirmed 
the intuitive view that there is a direct one-to-one correspondence between the degree 
of nonequilibrium present in the flow and the type of velocity distribution function 
needed to carry out the matching at a hybrid interface. 
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Because of the ’spiking’ of the pressure field in the two corners, one is not able 
to tell whether the well-known singularity in stress and heat flux that appears for 
the no-slip solution of the NS equations has been eliminated. Shown in Fig. 6 is the 
u-component of the velocity, the component that lies parallel to the lid itself. In this 
case, the lid is located on the near face and moves from left to right in the figure. 


u/c 0 



Fig. 6. A hybrid solution for the lid-driven cavity problem 
showing the component of the velocity lying parallel to the lid, 
for Kn = 0.01, Mud = 1.0 and a monatomic gas. The shearing 
lid is on the near face as indicated. 

The presence of velocity slip, seen in the figure as the difference between the plotted 
curve at the lid and the line at -1.0 indicating the Mach number of the lid, seems 
to have completely removed the infinite velocity gradient that otherwise appears in 
the two corners of the flow, and which leads to a singularity in stress when using the 
no-slip boundary condition. 

VI. Capability for 2D/3D Simulations 

Even though the lid-driven cavity problem represents two-dimensional flow, 
the interface used, consisting of a straight line, makes the interfacing effectively one 
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dimensional in terms of the spatial physics associated with the matching conditions. 
The next step in our work was to select a two-dimensional problem that required 
an interface that was also two dimensional. Such a case is shown in Fig. 7, which 
represents a supersonic flow past a rectangular prism. 



Fig. 7. A hybrid solution for Mach 16 flow over a rectangular 
prism, with NS shown in green and DSMC shown in blue. 


As this represents our first effort, the location of the interface was not positioned 
properly. The forward location of the interface should not place the NS solution 
(shown in green) ahead of the shock wave, but should confine it to the region closer 
to the boundary layer. However, the purpose of the test was to provide the experience 
needed to handle a more complex interface and to demonstrate that our overall pro- 
cedure was working properly. In this case the interface is a rectangular cutout with 
two sharp corners. The smooth transition seen in the figure between the NS solution 
and the DSMC simulation is what we had hoped to achieve before more complex 
testing was pursued. This work is being carried out by Duttweiler, as part of his 
thesis research, and a very preliminary report was presented by him at the CAS_96 
Workshop [7]. 

At an earlier stage of our work, Dahlby was able to test the upper limit of 
our capability to study true 3D flows using a pure DSMC simulation. He designed 
a generic SSTO vehicle body, which could be placed at an arbitrary angle of attack, 
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and ran a series of simulations on the Intel Delta parallel supercomputer. A total of 
1.1 million cells, 40 million particles, 512 nodes, 1000 sampling steps, and a run time 
of 110 minutes were used in each simulation. Figure 8 shows an example for 20° angle 
of attack, free stream Mach number of 10, and a monatomic gas. 



Fig. 8. Temperature distribution about a generic SSTO vehicle 
at 20° angle of attack, = 10, and for a monatomic gas. 

Less than half of the available memory was committed because of concern that poor 
load balance may cause a processor to be assigned a series of spatial blocks containing 
more particles than the processor had memory to store. This would cause the over- 
allocated processor to report a memory error and exit the simulation, which in turn 
would cause all other processors to abort the simulation as well. By using on average 
only half of the available memory, there would be a negligible chance of over-allocating 
a single processor. 

Dahlby was also able to conduct a number of interesting performance studies 
on our parallel code by using his generic SSTO body to represent a reasonably large 
problem. In one of his studies, he chose a simulation dimension of 120x120x120 for 
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a total of 1.7 million cells. This allowed him to evenly divide the simulation domain 
into cubical blocks of lengths 3, 4, 5, 6, 8, and 12 cells so that different runs could 
be carried out with a varying number of blocks per processor. Also, timing results 
were obtained with a total of 100 sampling time steps, as opposed to real runs of 
more than ten times greater. The aim was to study the effect on computation time 
of the tradeoff between the time required for load balancing, the final load balance 
state, and communication time. When there are only a few blocks per processor then 
the final load balance is poor because the apportionment of blocks is very restricted, 
and when the there are a large number of small blocks then communication cost 
and time to load balance become large. Figure 9 gives a summary of his findings. 
The most important observation made is that the histogram for the run-time data 
shows the minimum to be very broad, ranging from approximately 8 to 50 blocks 
per processor. The conclusion one draws from this is that, fortunately, most any 
reasonable engineering guess would turn out to be nearly optimum. 



Blocks per Processor 

2 4 6 8 10 12 14 

Block Length (cells) 

Fig. 9. Relative run time on the Intel Delta as a function of 
block size for the generic SSTO problem. 

Considering our current capability to carry out 3D simulations with a large 
number of cells and a correspondingly large number of particles, we were interested 
in testing the limit of our resolution capability for a typical 2D problem. A circular 
cylinder was chosen for this test using a wind tunnel dimension of 576 cells by 896 
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cells for a total of 0.52 million cells; the number density was 16 particles per cell on 
average for a total of 8.3 million particles. The free stream Mach number was set to 
10.95, the Knudsen number to 0.02, based on the cylinder radius, and the gas was 
monatomic. The temperature distribution about the cylinder is shown in Fig. 10. 



Fig. 10. Temperature distribution about a circular cylinder in 
a Mach 10.95 flow, monatomic gas and Kn = 0.02. 


This DSMC simulation was used by Dahlby in his thesis as a reference solution in his 
study of the effect of cell size on the position of the shock wave and the thickness of 
the thermal boundary layer on the cylinder. 


VII. DISCUSSION AND CONCLUSIONS 

As shown by the following list of publications, our progress during the reporting 
period has been considerable. Most of the goals outlined in the introductory discus- 
sion have been met and clearly suggest that further effort along these lines could be 
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productively pursued. On this basis, our ongoing work has dealt with (1) questions 
related to more complex geometries in two and three dimensions, (2) further testing 
of matching conditions for these more complex geometries, (3) steps that could be 
taken to further reduce the computational load represented by the sampling needed 
to create new particles at the hybrid interface, (4) the introduction of grid refinement 
in the Navier-Stokes portion of the hybrid, and (5) studies of computational time 
associated with different approaches so that the best method can be identified. 
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APPENDIX A 


Kinetic Flux- Vector Splitting for the Navier-Stokes Equations 

S. Y. Chou and D. Baganoff 
Department of Aeronautics and Astronautics 
Stanford University, Stanford, CA 94305 

(Submitted to: Journal of Computational Physics, Oct. 4, 1995; Revs: Feb. 26, May 24, 1996) 

Before a hybrid scheme can be developed combining the direct simulation Monte Carlo (DSMC) 
method and a Navier-Stokes (NS) representation, one must have access to compatible kinetic-split 
fluxes from the NS portion of the hybrid scheme. The kinetic theory basis is given for the devel- 
opment of the required fluxes from the Chapman-Enskog velocity distribution function for a simple 
gas; and these are then extended to a polyatomic gas by use of the Eucken approximation. The 
derived fluxes are then used to implement boundary conditions at solid surfaces that are based on 
concepts associated with kinetic theory and the DSMC method. This approach is shown to lead 
to temperature slip and velocity slip as a natural outcome of the new formulation, a requirement 
for use in the near-continuum regime where DSMC and NS must be joined. Several different flows, 
for which solid boundaries are not present, are computed using the derived fluxes together with a 
second-order finite- volume scheme and the results are shown to agree well with several established 
numerical schemes for the NS equations. 


I. INTRODUCTION 

A question that frequently arises when attempting to model plume flows generated by thrusters operating in 
the vacuum of space is: how to interface a numerical solution of the Navier-Stokes (NS) equations, which is most 
appropriate for the high-density gas flow found inside a nozzle, with the direct simulation Monte Carlo (DSMC) 
method [2], which is most appropriate for modeling the very low-density gas flow found in the outer plume? A similar 
question arises for the case of a blunt body in a rarefied high-enthalpy flow, where a numerical solution of the NS 
equations is often the most appropriate for handling the high-density gas layer found near the relatively cold body 
surface, while the DSMC method is the most appropriate for the rarefied outer flow. 

Both questions can best be understood by considering the much simpler case of an impulsively started piston in a 
stationary gas, where the piston travels in a direction normal to its surface, the piston Mach number is assumed high, 
and the piston temperature and gas temperature are initially equal. A NS solution for a representative case is shown 
in Fig. 1, where the piston Mach number is taken to be 5.0 and the gas to be ideal and monatomic. A normal shock 
wave forms ahead of the piston, the gas temperature (dashed curve) rises behind the shock and then falls as a result 
of the cold piston surface (located at x/L — 1). Because of the nearly constant pressure between the shock and the 
piston, the density (solid curve) rises in the thermal layer in an inverse ratio to the temperature. In this case, the 
peak density is 8.2 times the density behind the shock wave. Although the DSMC method must be used to obtain 
the proper shock wave profile, as is well known, it is clear that the much higher density in the thermal layer would 
lead to a greatly increased DSMC simulation cost in that region, prompting consideration of a hybrid scheme, where 
the DSMC method would be used to model the outer flow and the NS equations to model the thermal layer. This 
concept is schematically depicted in Fig. 2. 

Because the DSMC method is based on kinetic theory, the DSMC fluxes to be matched at an interface of a hybrid 
solution are physical quantities: the fluxes emanating from the DSMC region F + are one-sided kinetic fluxes and the 
frame of reference is an inertial frame, as indicated in Fig. 2. The flux of mass (momentum or energy) from the DSMC 
side is gotten directly from the particles that cross the interface in one time step. In order to develop an effective 
hybrid scheme, one must therefore address the following issues. 

• The same definitions should be used for the split fluxes when matching at an interface, preferably the kinetic 
theory definition of one-sided kinetic fluxes (DSMC) in an inertial frame of reference. 

• As indicated by Figs. 1 and 2, considerable nonequilibrium may exist at an interface in the torm of heat flux 
or shearing stress when considering a boundary-layer flow. Therefore, the numerical scheme chosen to interface 
with DSMC must solve the NS equations; and because the DSMC method possesses characteristics of a time- 
dependent finite- volume scheme, compatibility suggests the use of the same scheme for the NS portion. 

• Because matching can only be carried out in the near-continuum regime, where both DSMC and NS are valid, it 
is clear that velocity slip and/or temperature slip would be present at the body surface, i.e., the piston surface 
in Fig. 2. Therefore, the NS solution must account for slip; the no-slip boundary condition cannot be used. 
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• The numerical scheme handling the NS portion must be shown to agree with existing numerical methods for 
the NS equations before a hybrid scheme is constructed. 

• Conversion of the NS values for the one-sided kinetic fluxes at an interface into a corresponding collection of par- 
ticles for insertion into the DSMC simulation must be carried out in a physically compatible and computationally 
cost effective way. 

• Conditions must be identified that define a suitable location for placement of an interface. For example, for 
moderate nonequilibrium, it is generally found that good results are obtained for DSMC simulations if the local 
cell Knudsen number (local mean free path length to cell dimension) is greater than unity. To use this criterion, 
questions of compatibility with the NS equations would have to be explored. 

• Both the NS scheme selected and the resultant hybrid scheme created must be compared with different DSMC 
solutions for complete verification. 

The objective of the present study is to address and discuss the first four issues listed above; the remaining three will 
be covered in follow-on reports. 

Much of the groundwork for the present investigation can be found in a book by Patterson [13] and in papers 
by Broadwell [4], Bird [2], Pullin [15], Prendergast & Xu [14], Mandal & Deshpande [12], and Gooch [8]. Because 
we seek a formulation based on the NS equations, and because much of the previous work focused on the Euler 
equations, it is appropriate to first turn to general results found in work on kinetic theory, in particular, the well- 
known system of moment equations formed from the Boltzmann equation. This approach treats the equilibrium and 
nonequilibrium components of the fluxes on the same footing, and the derivation itself provides the definitions along 
with the justification for their use. An alternative is the separate introduction of the nonequilibrium terms at the 
macroscopic level, as discussed by Macrossan & Oliver [10] and by Mallett et al. [11], but this approach would not 
directly satisfy the first and fifth condition in the above list of issues. 


II. MOMENTS OF THE BOLTZMANN EQUATION 


Many of the equations and concepts to be presented in this section can be found in standard work on kinetic theory, 
such as Chapman & Cowling [5], Grad [9], Patterson [13], Vincenti & Kruger [20], Woods [22], and Bird [3]. We start 
with the case of an ideal monatomic gas in the absence of external forces and assume the gas is sufficiently dilute for 
binary collisions to dominate. For this case, the Boltzmann equations reads 


d{nf) d(nf) _ 
dt + Ck dx k 


d{rif ) 


dt 


( 1 ) 


zoll 


where n is the number density, / is the velocity distribution function, c k the molecular velocity in an inertial frame, 
the repeated index k denotes a sum, and the right-hand side represents the collision integral. The moment equations 
are obtained by multiplying the Boltzmann equation by any function of molecular velocity Q(c*) and integrating over 
velocity space. These equations are represented by 


^ (9 

^-(n < Q >) + (n < c k Q >) — A[Q]. 
The two operators appearing in (2) are defined by 

/ OO rOO pOO 

/ / Qfdcidc-zdcz 

- OO J —OO J — OO 


and 


/ oo poo poo r 

/ Q \ 

-oo J — oo J — oo L 


d(nf) 


dt 


dcidc2dc3 . 
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( 3 ) 
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When the arbitrary function of molecular velocity Q(c,-) is chosen to be one of the five collisional invariants Q = 
m {l, c,-, c 2 /2}, where m is the molecular mass and c 2 represents the square of the velocity magnitude, then the 

/ pj y 

corresponding moment of the collision integral is identically zero, i.e., A[Q ] = 0. This is a general result that 
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holds for any distribution function / and for any molecular interaction law. This selection leads to the conservation 
laws for gas dynamics, which can be written in the form 


d . ,'NV d . „ INV 

— (n < Q >) + ^-(u < CfcQ >) = 0, 

(5) 

or, when using each of the collisional invariants in turn, one obtains the set 


d d 

Wt (p) + d^ {p < Ck >] = 0 

(6a) 

-Q t i P < a >) + ^( p( c k c i » = 0 

(6b) 

|(p < cV 2 » + £(, < c ,cV2 » = 0, 

(6c) 

where p — rnn is the mass density. 

Introduction of the thermal velocity components C; = (c, — u*), where u, =< c,- 
allows one to introduce the central moments defined by 

> is the mean or fluid velocity, 

Pij — p < C{Cj > 

(7) 

V - Pkk/% 

(8) 

T{j — —Pij + P$ij 

(9) 

e =< C 2 / 2 > 

(10) 

9i = p < CiC 2 / 2 >, 

(11) 


where Pij is the pressure tensor or stress tensor, p the pressure, is the viscous stress tensor, e is the internal energy 
(translational) for a monatomic gas, and <?,• is the heat flux vector for a monatomic gas. The conservation laws for 
gas dynamics can then be written in the familiar form 

§t {p) + £; {pUk) = 0 (12a) 


ft ft 

dt ^ pUi ) + (P UkUi + Pk '} ~ 0 


(12b) 


d 

f “ 2 M 

d 

( u2 \ 

dt 

K £+ t)J 

+ dx k 

puk 1 e + ~2 ) + Pk ' Ui + qk 


(12c) 


If the gas is not simply a monatomic gas but has internal structure, then the above procedure must be modified. 
Because the general problem includes the question of what equation replaces (1), the problem is rather difficult and 
it becomes necessary to make use of a suitable approximation. One simplistic approach is to assume that all internal 
molecular energy modes are in equilibrium, both internally and with the translational degrees of freedom. Thus, 
the additional internal energy ej nt can be expressed in terms of the translational temperature T by the equilibrium 
relation 


_ 1 
e int - ^ 


5 - 3 7 
7 ~ 1 


RT, 


( 13 ) 
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where R is the gas constant, and where an accounting for the additional internal energy is introduced through the 
ratio of specific heats 7 , and for which e, nt = 0 in the case of a monatomic gas. 

To properly account for the amount of energy that is carried by a particle with internal structure, the energy me 2 / 2 
must be replaced by (me 2 f 2 + e), where e is the additional internal energy per particle, and therefore, the quantities 
of interest become 

Q = {m, me,, (mc 2 /2 + f)} . (14) 

Assuming Eq. (1) continues to hold for the extended distribution function /(c,-,e) and provided the right-hand side 
is interpreted in a suitable way, then an additional integral over e is formally required in applying both (3) and 
(4). It is reasonable to argue, however, that the quantities in (14) must continue to be conserved in a collision, and 
consequently, (4) again evaluates to zero; and thus, (5) remains unchanged. 

In evaluating the left-hand side of (5) for the five different quantities in (14), identical results to those obtained for 
the monatomic gas will be found for all quantities that contain polynomials in c» alone. This follows from the fact 
that integration over the e variable can be taken first and independently from the <7 integration. Therefore, Eqs. ( 6 a) 
and ( 6 b), and consequently (12a) and (12b), are fully recovered. The same conclusion also applies to the first term in 
the quantity (me 2 / 2 + e) and therefore ( 6 c) is replaced by 

— (p < c 2 /2 > +n < c >) + -z — (p < c k c 2 / 2 > +n < c k c >) = 0. (15) 

ot OX k 

Substitution of the central moments (7)-(ll) into (15) reproduces all of the previous algebra for the monatomic gas 
and leads to 

d ( u 2 \ d ( u 2 \ 

— p ( e + — J - 1 - peim + 7 ^- puk \e + — J + PkiUi + q k + (n < C k £ > +pukei n t) = 0 , (16) 

where the definition < e >— mei nt has been used. Therefore (12c) also continues to hold provided we replace (10) by 

e = (< C 2 /2 > +e int ) (17) 

and the definition of the heat flux vector ( 11 ) by 

qi = p < CiC 2 / 2 > +n < Ci€ > . (18) 

We therefore conclude that the complete collection of equations (5), ( 6 ) and (12) can be used as they stand, provided 
definitions (17) and (18) are employed when the gas possesses internal structure, and a state of equilibrium exists 
between the internal modes and the translational degrees of freedom (see [20], p. 326). 

Because the conservation equations (12) can be developed for any general fluid through use of phenomenological 
arguments alone, the set is actually more general than the kinetic theory derivation would indicate, i.e., they are also 
the conservation equations for fluid dynamics. However, we are only interested in treating an ideal gas flow and the 
use of the kinetic theory approach is necessary because it shows that the set is valid for any degree of translational 
nonequilibrium, that is, for any translational velocity distribution function one cares to consider. If one chooses the 

Max 

equilibrium distribution, namely the Maxwellian distribution / , then the set becomes the Euler equations, because 

viscous stress and heat flux are identically zero for / . This step allows one, in effect, to reproduce the special 

case considered by Pullin [15] in his equilibrium flux method (EFM), as well as the work by Mandal & Deshpande 
[12], on kinetic flux- vector splitting (KFVS) for the Euler equations. The inviscid limit of the recent work by Mallett 
et al. [11] is also recovered by this step. If one chooses a Chapman- Enskog (CE) distribution / , then the set 

becomes the Navier-Stokes equations, because stress and heat flux are then given by the corresponding Chapman- 
Enskog expressions. This is the path that we will follow in developing a KFVS scheme for the NS equations. One 
may also choose a discrete representation for } and the equations still hold. This is an alternate interpretation of the 
concept that lies behind the state vector splitting scheme for the Navier-Stokes equations introduced by Gooch [ 8 ]. 
The most important point is that one is free to choose any translational velocity distribution function whatsoever in 
using Eqs. (5), ( 6 ) or (12); and in so doing, the set becomes closed, as long as / is fully specified. Otherwise, if j 
remains general, then one is faced with the well-known closure problem, when using a moment method, because r,y 
and qi are then unknown quantities in the equations. It is useful to emphasize the fact that the conservation equations 
as displayed in (12) are not the NS equations until one introduces / 

Each of the five separate moment equations ‘represented by either (5), ( 6 ) or (12) can be expressed by the single 
form 
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(19) 


8U_ dF\_ 
dt dxk 


= 0. 


The complete specification of Fi in three dimensions, as defined by (19), requires the evaluation of 15 quantities. But 
the task is made considerably simpler, if one considers a finite volume scheme, uses Gauss’ divergence theorem, and 
writes Eq. (19) as 


If 

dt J v 


UdV + [ F n dS = 0, 
v Js 


( 20 ) 


where S encloses the volume V and F n is the projection of Fi onto the unit outward pointing normal for the surface 
element dS. If V is taken to be a rectangular volume, then one only needs to evaluate five quantities for each planar 
surface, a conceptually simpler task provided F n can be evaluated directly. Using the notation of (5), we then have 


INV 

U = n<Q > 


( 21 ) 


and 

l N V 

F n = n < c n Q >, 


( 22 ) 


where U is the state vector, F n is the total flux vector, and c n is the component of the molecular velocity normal to 
the planar surface. Equation (22) clearly represents a physical concept. This can be seen from the fact that Q is a 
scalar quantity which is carried across the fixed surface by c n , thus creating a physical flux in that quantity. 

The five fluxes defined by (22) are total fluxes, not the one-sided fluxes depicted in Fig. 2. In addition, these general 
expressions contain both the inviscid fluxes as well a s the nonequilibrium components due to viscous stress and heat 
flux, as can be seen from the corresponding terms in (12). Because many quantities evaluate to zero in arriving at the 
set (12), and these cannot be recovered in a simple way, one must use the more primitive set (5), or (6), in developing 
the algebra for the one-sided, kinetic-split fluxes. 


III. KINETIC SPLIT FLUXES 


As seen in Fig. 2, we are interested in expressions for the one-sided fluxes based on a fixed interface and an inertial 
frame of reference. On this basis and on considering the xi direction as positive, we can split the integration in ci as 
follows 


/ oo poo p oo p oo poo / pO pea \ 

/ / {....} dcidc 2 dcs = / / ( / + / \ {....} dc\dc 2 dc 3 

„oo J — OO «/ — OO «/ — OO J ” OO — OO J 0 / 

/ OO pOO / /* — til f 00 \ 

/ / +/ {....} dC x dC 2 dC 3y 

- OO J — OO — OO j —U 1/ 


(23) 


where the second expression introduces the thermal velocity components C, ■ The concept of a kinetic split flux has a 
long history in kinetic theory and a clear application to an equilibrium flow can be found in a textbook by Patterson 
(Ref. [13], pp. 163-167). On using the following notation to represent the splitting in (23) 


F = F~ + F + , 


(24) 


and on introducing a Cartesian coordinate system (n,tl,t2) located in an arbitrary fixed planar surface, we obtain 
from (22), (23) and (24) the definitions 

/ OO r OO /> — ti n 

/ / (C„+u n )Q'" v fdC„dC n dC t2 (25) 

-OO J — oo J — oo 


/ oo poo poo 

/ / (G, + «„)(? fdC„dCndC t2 (26) 

-oo */ — OO j — ti n 

These relations provide the means for computing F n when / is a known function. 

Because the total fluxes are known from the terms in (12), explicit expressions for each of the collisional invariants 
given by (14) need only be listed for, say, F+. It is also useful to introduce a more physically descriptive notation for 
the split fluxes as follows. 
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fdC n dCtidC t2 


(27) 


F + 

A zero 



F+ 

mass 



( C n + u n ) f dC n dCtidCt 2 


(28) 


F + 

n — mom 



( Cn + u n ) 2 fdC n dCndC t2 


(29) 


/ OO poo p oo 

/ / (C„ + u n )(Cn + un)fdC n dCndC t2 (30) 

-OO J — oo J — u n 

/ oo poo poo 

/ / (Cn + «„) [(Cn + «n) 2 + {C t l + U t \) 2 + (C* 2 + U t2 ) 2 ] 

-oo J — oo J — u n 

x l -fdCndCnd.Cn (31) 

p OO pOO poo poo 

F,t,-er, ergy = n / / / (C n + u n )cfdC n dC tl dC t2 de = A,+„ cie „ + />u n e+, (32) 

JO J — oo J — oo J — u„ 


p+ — zr+ 4- P.+ 

- 1 energy x tr— energy 1 int — energy' 


(33) 


Integration over the e variable is not shown in Eqs. (27)-(31) because it can be carried out as an independent operation 
and done first. However, it does appear explicitly in (32), and this equation must be included when a gas has internal 
structure. The role played by (32) can be seen from the fact that the C n e term in the integrand is the source of the 
second term in (18), while the u n e term is the source of the second term in (17). The difficulty in applying (32) results 
from the fact that we need an explicit expression for the joint distribution /(C,-,e) in order to evaluate these split 
fluxes; and this lack of knowledge is indicated by the notation employed following the second equality. However, in 
the case of the total e int , it is given by (13) since our model assume equilibrium for the internal degrees of freedom. 
Likewise, for the total flux represented by Ag Bucfcen , we can use the Eucken model which replaces < C n t > by a 
quantity that is proportional to the temperature gradient, thus making it proportional to the heat flux vector (see 
[22], p. 66]). The net effect changes the value of the coefficient of thermal conductivity, as well as the Prandtl number, 
from that for a simple gas to the proper values for a gas with internal structure. This then allows one to use the 
same basic relations found for a simple gas. We will use the Eucken approximation and represent the incremental 
contribution to the heat flux due to the internal degrees of freedom by 




= 




= -K VT. 


(34) 


These steps alone, however, still do not answer the question of how one evaluates the split fluxes appearing in (32). 
This will be done after we have completed the implementation of the Eucken model. 

The extra quantity F+ ro is also listed because it defines how the velocity distribution function itself is split, which 
is needed in the overall algebra. This can be seen from the normalization condition F zero — F~ r0 + F+ er0 = 1 for a 
probability distribution. The expression for F+ 2 -. mom is not listed because it can be gotten by merely interchanging 
fl and t2 in (30). For a monatomic gas, the set of equations is very general and applies for any velocity distribution 
function one cares to define, for example, even for a discrete distribution. For a polyatomic gas, the set is not quite 
as general and is limited by the Eucken model and the approximations to be introduced below. Our interest is in the 
Chapmann-Enskog distribution and the resulting split fluxes. 


IV. CHAPMAN-ENSKOG SPLIT FLUXES 

A gas flow that is in thermodynamic equilibrium is represented locally by a Maxwellian distribution, and a gas flow 
that is slightly disturbed from the equilibrium state is represented locally by the Chapman-Enskog distribution. The 


6 



CE distribution is obtained as an approximate solution of the Boltzmann equation (for a simple gas) and is expressed 
as a product of a local Maxwellian and a polynomial function of the thermal velocity components C,, that is, by the 
relation 

/ —f (l + ^i + ^ 2 ) (35) 


where 


f MaX = (2TrRT)~ 3l2 exp(-C 2 /2RT) 

* = -(?) 


and where K (1) is the coefficient of thermal conductivity and p ) the coefficient of viscosity as determined by the 
first-order Chapmann-Enskog procedure, and Sj k the Kronecker delta. Because both the temperature gradient and 

the velocity-gradient tensor appear as parameters in / , notational efficiency can be gained in the algebra that 

follows by replacing those quantities by the Chapman- Enskog expressions for stress and heat flux, i.e., by 


Qi 


CE 


= -K 


(i) 


&T 

dxi 


(36) 


ce (i ) ( dui duj\ 2 (''>{du k \ 

T « = " + (37) 

At this point we face a logical difficulty. Must we fix the Prandtl number to the value for a simple gas (Pr = 
/i (I> c p / K W = 2/3) or may we allow it to vary so that (36) becomes consistent with the Eucken approximation? 
As there is no simple answer, we will defer the question to the point in the analysis where logical conflicts can be 
more easily identified. On substituting (35) into Eqs. (27)-(31) and using the same orthogonal coordinate system 
(n,tl ) t2), we see that in each case the integrand becomes a product of polynomials in the thermal velocities Ci and 
the Maxwellian distribution f M “ x . The only difficulty that appears in carrying out the integration is the large number 
of terms that are produced. For example, <j)\ and 0 2 are composed of 4 and 9 terms, respectively, and therefore (35) 
leads to a total of 14 terms. On evaluating F^_ energy alone, we are faced with 252 terms. This at first appears to 
be an overwhelming task, until it is noticed that many of the terms are zero, because in the C t 1 and C t 2 variables all 
odd moments of the symmetric function / are zero and even moments are well-known functions of RT. Likewise, 
integration in the C n component can be handled by splitting the integration into — ttoo to 0 + 0 to 00 which then 
leads to exponential and error functions (see [3], p. 417). Because the collection of functions is small, there is hope 
the results can be finally assembled into fairly compact expressions. Even though the details are daunting, Patterson 
(Ref. [13], p. 77) used a method he reports was ’’initiated by Maxwell and developed by Chapman” to develop 
general expressions for the total fluxes, for the case of nonisentropic flow, which employs concepts and algebra very 
similar to that employed here. Our interest, however, is in obtaining the split fluxes for which many more terms 
must be handled. To carry out the present study, intensive use of symbolic mathematical manipulations, provided 
by MATHEMATIC A, was made in handling the algebra. Only the final collected integrated results will be presented 
here. A discussion of the detailed steps required to obtain the following relations can be found in [6]. 

FLo = \ [(1 ± «i) ± «2 (S„t™ + (2Sl - l)C E )] (38) 

FL„ =Px/ra72[(l±ai)S„±a 2 (l-xi)] (39) 
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± «2 ( S n + q„ E ) 


(40) 
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V^T [S„F±„„] + \p [-(1 ± ± a 2 qg E ] 


(41) 
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tr — energy 


ps/mfi (l±oi) (s„(^+S 2 ) + x 2 ) ±a 2 (2 + S 2 +* 3 ) 


(42) 
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r int — energy 


( a ?L fcen +/ ?u « e j fc «t) 



(43) 
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(44) 
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X2 = 
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erf(S n ), 
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+ -f CE 

T 2 >nn 
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+ S t 2Q t 2 ” Xi (1 + S tl + S t2 ) - 

u n /V2RT , 


- CE /p 

nn l r) 


S' ! = si + s, 2 , + S'i 2 

2 


= ^?? E /(pV5fir). 



Each individual component of 5, , fjj and qi is not listed, as they are clearly nondimensionalized the same way. 

Equations (38)-(42) are not outwardly affected by the physics associated with additional internal degrees of freedom, 
i.e., 7 ^ 5/3. In view of this, these equations should not contain 7 explicitly, when expressed appropriately. This is 
easily done by introducing the speed ratio S — u/y/2RT, which is frequently used in kinetic theory, as opposed to 
the Mach number, which requires the introduction of 7 through the isentropic speed of sound. Use of the speed ratio 
S not only provides a useful physical check on the mathematical results, it also allows the final expressions to be 
written in a more compact form. On the other hand, the physical concepts that lead to (43) directly involve additional 
internal degrees of freedom and the relation should contain 7 explicitly, which is seen in the expression that follows 
the second equality, a step to be explained in the following discussion. 

It is appropriate to review several consistency checks on Eqs. (38)-(43). First, the total fluxes (24) should agree 
with the corresponding expressions in (12). Starting with (38), we have F zero = 1, which is the correct normalization 
condition for a probability distribution. From (39)-(43), we have 


Fmass — P^n 


(45) 


Fn-, 


— ™,2 


P U n + P ~ T n 


(46) 


Ft\— mom — P^n^tl ^ntl 


(47) 
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(48) 


Fi n t-energy — {^Qsuck^n + pUn e int) 


1 

2 


5-37 
7 - 1 


PU n - 


(49) 


These correspond to all the expressions in (12) for the orthogonal coordinated system used. Equation (49) provides 
both the extra term required by (17) as well as the extra term required by (18). However, this returns us to the 
logical difficulty raised in the discussion following (36) and (37). If the Eucken approximation is introduced when first 
using (36), then its effect will appear twice when summing (48) and (49), i.e., in the combination ( q n + A q Euckcn )- 
Knowing this, one approach would be to introduce the Eucken approximation here by absorbing A q Euckcn into the 
term q n by using Eq. (34). Although it appears to be a reasonable step, it is in fact a bold step, because split fluxes 
are required in (42) and (43) and one cannot be sure that when the Eucken approximation is introduced into Xj and 
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X 3 in (42) that it will properly account for the absorption of the split quantities A q^ uckcn in (43). Consequently, one 
may encounter incorrect energy split fluxes at an interface or a boundary. Actually, we have no alternative, as we do 
not have /(C,-, e) with which to compute these split quantities. In the following development, we will use the Eucken 
approximation throughout and assume that proper accounting is made for A</* ucfcen by and X 3 - 

An additional assumption is actually needed to complete the specification of (43), and this involves the evaluation 
of ef nt , which again requires knowledge of /(C,-, e). Assuming the Eucken model properly accounts for the correlation 
Ag s ucfcen = n < C n (. > then it may be permissible to employ the equilibrium assumption to approximate ef nt . The 
assumption that the internal energy modes are in equilibrium, both internally and with the translational degrees of 
freedom, leads to the conclusion that Ci and e are statistically independent random variables, and therefore, /(C t , e) 
reduces to a product function. On this basis ef nt can be evaluated, and this step leads to the expression following 
the second equality in (43). However, we should note that if the same assumption were used to evaluate < C n e > as 
well, then we would have A q Eucken = n < C n e >= n < C n >< <f >= 0 because < C n >= 0 by definition; and this 
would lead to the loss of the Eucken approximation. In summary, on combining (42) and (43), our approach consists 
of absorbing ^q^ uckcn into \2 and X3 by the introduction of the Eucken approximation and using the equilibrium 
assumption to evaluate ef nt ; this leads to the second expression in (43) as well as the second expression in (49). 

A second check leads to the requirement that one should recover the known values for the Maxwellian distribution 


when setting the nonequilibrium parameters f and q to zero. 

= 1(1 ± ai ) (50) 

FL„ =Pv/H772[(l±a,)S„ ± 02 ] (51) 

= P [(1 ± «i) (^ + i) ± a 2 S„] (52) 

Fti-mon, = V2RT *.„] (53) 

(l±a,)S„(5 + S 2 )±o 2 (2 + S 2 ) (54) 

FL-^ rgy = \ (~y) [RTF^„] . (55) 


These expressions are in full agreement with results obtained by Patterson (Ref. [13], Section 5.3, equations (4), (11), 
(17) and (19)) where his objective was to compute the momentum and energy exchange at a surface for an equilibrium 
gas. These relations were later introduced in work by Pullin [15] (for one-dimensional flow) and more recently by 
Mandal h Deshpande [12] and by Mallett et al. [11]. 

Finally, on setting the fluid velocities to zero ( S — 0,ai = 0,<*2 = l/V^F), one must recover the well-known kinetic 
theory values for a stationary equilibrium distribution (Maxwellian) 


Ft„c = 1/2 

(56) 

/mass ~ ±P\JRT/‘2 k 

(57) 

— morn iP/2 

(58) 

^U-mom ~ 0 

(59) 

F^-energy = ±P\/2^T/7r 

(60) 
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(61) 


PL-, „„„ = ±5 (^=t) ps/m^- 

The first five relations clearly agree with known values in kinetic theory, while the final relation is based on several 
critical assumptions as discussed above. 

The importance of the nonequilibrium parameters f^ k and q k in Eqs. (39)-(44) can be judged by comparison with 
equivalent Euler split fluxes. Of interest in a comparison are the Steger-Warming split fluxes [17] and the equilibrium 
values given by (51)-(55). Obviously, the nonequilibrium parameters are not constants in a flow, but in a graphical 
display representative values are quite useful, and a representative peak value for both parameters in a moderately 
strong normal shock wave is roughly -0.3. This value along with 7 = 1.4 and the assumption of a one-dimensional 
flow were used to develop Fig. 3, which presents the momentum and energy split fluxes versus the local Mach number 
for three cases: KFVS for the NS equations given by (40) and (44); KFVS for the Euler equations given by (52) and 
the sum of (54) and (55); and Steger-Warming splitting. On the scale shown, the results for the equilibrium values 
and Steger-Warming group fairly closely together (a comparison apparently first made by Mandal and Deshpande 
[12]), while the present work shows considerable difference, especially in the asymmetrical shift seen in F + and F 
The figure clearly point to the fact that in a hybrid solution, where the one-sided fluxes are to be matched at an 
interface having nonequilibrium conditions, the present kinetic split fluxes must be used. This is because one does not 
otherwise know how to adjust the Euler split fluxes, to account for viscous and heat conduction effects, even though 
the total fluxes are known from ( 12 ). 


V. BOUNDARY CONDITIONS 


When a nonreacting particle in the DSMC method passes through a body surface during a time step, the procedure 
is to emit the same particle from the surface with a new velocity depending on the boundary conditions. Thus, the 
DSMC method effectively treats a solid boundary as though it too consists of a gas, but at different conditions. Ihis 
concept was clearly described by Patterson (Ref. [13], p. 165) and used in his analysis of molecular interactions with 
boundaries. Because the same particle is emitted, the wall gas is identically the same gas, and therefore the two share 
the same molecular mass and gas constant. Also, because every particle passing through a body surface is treated 
this way, the number of particles per unit time per unit area passing into a wall is exactly balanced by the rate of 
emission. However, this condition does not tell us what the number density of the wall gas is, i.e., it does not fix n w , 
or p W) for the wall gas, nor does it fix the temperature of the wall gas T w . What we do know is that the total mass 
flux must be zero at a material surface, i.e. 


(F m ass) surf ace — (^majj)g 4" — (®2) 

where gas and wall positions, consistent with the end-wall geometry shown in Fig. 2, are assumed. Generally speaking, 
in the DSMC method one often assumes that the wall gas is in equilibrium and the wall is stationary. On this basis, 
Eq. (57) can be used in (62) to describe the wall gas, and we can write 

= ( 63 ) 

A particularly simple boundary condition is the case where the wall temperature T w is specified, namely, an 
isothermal boundary condition; and for this case, Eq. (63) fixes the density of the wall gas p w , since (F+ a3S ) g would 
be known from (39) and the state of the gas from the update procedure for the numerical solution of the NS equations. 
Given p w , T w and p w = p w RT Wy we then have on using (58) the relation 

(F n — mom) surf ace — (^ + - mom ) 3 +P„/2, (64) 

which fixes the normal stress at the surface, since (Fn_ mom ) g likewise represents a known quantity from (40) and the 
update procedure. Likewise, on using (59), we have 


(Ft 


tl-mom ) surf ace 


(F'tl — mom) 9 ’ 


(65) 


which determines the tangential stress on the surface. Finally, the energy flux to the surfaces is found by using (60) 
and (61) together with (44) and the state of the gas from the update procedure, i.e. 


(F< 


energy) surf ace — (F^ n ergy)g P w \/2 RT W j 7T 




5 — 37 


( 66 ) 
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which completes the specification of the conditions at a surface for an isothermal wall. The most important outcome 
from this analysis is that the temperature of the gas near the surface may not be equal to the specified wall temperature 
T w , which gives rise to the possibility of temperature slip. Likewise, the tangential velocity near the surface may not 
be zero, leading to velocity slip; however, the normal velocity near the surface must be zero because the expressions 
in (62) also define the fluid velocities for the two separate gases. 

For an adiabatic boundary condition, we know that the total energy flux to the surface is zero, while the wall gas 
temperature and density are unknown. Thus, we have 


{F energy) sur / ace — (-^energy)? “1” (F e nergy)™ — 

and on using (60) and (61) for the wall gas, we have 

(■^energy)g = Pw \/2 RT W / K 1 + — ^ ^ 


(67) 


( 68 ) 


which along with (63) provides two equations in the two unknowns p w and T w , since {F+ ass ) g and (F+ nergy ) g are 
both known from the update procedure at each time step. Likewise, the normal stress on the surface can be obtained 
using (58) and the tangential stress using (59), completing the analysis. 

The application of flux boundary conditions to the wide variety of possible boundary conditions is quite straight- 
forward and only two are given here. For example, accommodation coefficients for momentum and energy are often 
used in the application of the DSMC method, and it would be a simple matter to include them as well. In addition, 
analytical results for both temperature slip and velocity slip can be obtained from Eqs. (62)- (68) and these can be 
shown to reproduce the special case considered by Patterson (Ref. [13], Section 4.4, equations (30)-(33)). A detailed 
discussion of this topic will be covered in a follow-on report. 


VI. NUMERICAL COMPARISONS 


Our main objective is to show that the derived expressions for the split kinetic fluxes given by Eqs. (39)- (44), together 
with the flux boundary conditions, lead to valid solutions of the NS equations. For this purpose, comparisons are 
made with first-order schemes by Steger k Warming [17] and Roe [16] , while second-order comparisons are made with 
the symmetric limited positive second-order (SLIP2) scheme reported by Tatsumi, Martinelli & Jameson [18]. Central 
differencing was used in each of these cases to handle the terms introduced by viscous stress and heat flux. A version 
of the monotone upstream-centered scheme for conservation laws (MUSCL), van Leer [19], was used in applying the 
KFVS relations for the NS equations. The parameters chosen in using MUSCL are defined in the following relations 


t?MU SCL r '+ i T?~ 

F , + i - F + i + 


where 


C +i = F + (U!- +i ) 

F~ +i = F-(U« +i ) 

n 


i + 
7~ 

‘+i 

IJ L 


AU- + , = U i+1 
minmod(a, b ) 


A Ui_ h ^L 

*6 r 

i = U ,+ 1 — —minmod A f/ <+ |, 


Ui + -A minmod 



M < 1*1 
M > 1*1 


(69) 


and where the L/R arrangement is defined in Fig. 4. A first-order scheme is obtained for 9\ — 0 and a second-order 
scheme for 6\ = 1. Most of the work was carried out for f3 — 1.5. The function F+(U L ), for example, represents any 
one of the equations (39), (40), (41) or (44). These equations contain spatial derivatives of U as well, and the concept 
of the function must be generalized to include them. For a problem with one spatial dimension, the update formula 
(first-order in time) is then given by 

yn + l = {J n _ f MUSCL _ F MUSCL\ n > ( 70 ) 

’ ’ AX \ ,+ 2 * 2 ) 
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where n represents the time step. 

The problem of determining the profile for a normal shock wave represents a steady flow for which viscous stress 
and heat conduction are very important; and this represents an extreme limit that provides appropriate conditions for 
testing. Although it is well known that the NS equations do not give physically realistic profile predictions for strong 
shock waves, we are only interested in establishing that valid NS solutions are being obtained. For this purpose, we 
have the method suggested by von Mises [21] and by Gilbarg &; Paolucci [7] for solving the NS equations for steady 
one-dimensional flow, which can be used as a standard of comparison. To obtain a reliable reference, a four-step 
Runge-Kutta method and 1200 points in the domain of integration were used. The primary integration takes place 
in a computational space where the fluid velocity is the independent variable and, because the normal stress and the 
heat flux variables seem to show the greatest numerical sensitivity, these were selected for display. Figures 5 and 6 
give results for a shock-wave Mach number of 1,5 and air as a representative gas. First-order schemes for KFVS, 
Steger-Warming and Roe are compared with the reference solution (solid curve) in Fig. 5; and it is clear that Roe’s 
scheme compares the best, with Steger-Warming exhibiting the familiar transition near the sonic station and KFVS 
showing only moderate success. Second-order schemes consisting of KFVS and Jameson’s SLIP2, along with the 
reference solution, are compared in Fig. 6; and it is clear that it is not possible to distinguish between them. Higher 
Mach numbers were also studied and the same general observations were made. 

In the same spirit of reviewing extreme conditions, it is appropriate to consider a flow that corresponds to the 
Euler limit, namely, a shock-tube flow. Because of the presence of the contact surface, it is most useful to display 
the density and temperature variables as they frequently exhibit large changes there. Figs. 7 and 8 present results for 
pressure ratios (driver to driven gas) of 3 and 20, respectively, and where the flow is from left to right. Starting from 
the left, the transitions seen are the expansion fan, the contact surface, and the shock wave. We see that there is good 
agreement between second-order KFVS and SLIP2, with the largest difference at the contact surface, possibly a result 
of the different limiters used. Roe’s first-order scheme clearly shows severe rounding of the profiles. The separation 
distance between the shock wave and the contact surface grows linearly with time and therefore the relative rounding 
of the corners would appear to diminish with increasing time. An early time was specifically chosen for display to 
emphasize the relative difference between the different schemes. 

In both the ideal shock-tube and shock-wave problems solid boundaries are absent. However, the case ol an 
impulsively started piston is an example where nonequilibrium effects near a solid boundary can lead to temperature 
slip. Figure 9 shows an expanded view of the thermal layer near an isothermal piston, where the piston Mach 
number is taken to be unity and the gas is ideal and monatomic, and it is clear that the gas temperature near 
the piston is higher than the piston temperature, as a result of the flux type (DSMC) boundary conditions used in 
solving the problem. A theoretical expression for temperature slip in a slightly rarefied flow of a monatomic gas 
was developed by Patterson [13] (see equation (33), p. 125) and his prediction is shown by the circle symbol. The 
excellent agreement seen undoubtedly results from the fact that Petterson’s theory makes use of the assumptions in 
NS, and therefore, boundary conditions (62)-(68) used in the numerical solution and the analytical approach used by 
Patterson correspond closely. A necessary next step, of course, is to carry out comparisons with DSMC simulations 
to determine the conditions under which the magnitude of the jump is physically correct. This study, which requires 
careful discussion, will be presented in a follow-on report. The most important observation is that the method 
introduces slip into the NS formulation in a very natural way through the use of flux boundary conditions. 

A final problem to be reviewed is the steady, two-dimensional flow produced by a plate sliding across the open end 
of a square cavity, and for which the temperatures of all four material surfaces are held fixed and equal. This is a 
well-known problem for which the NS solution, based on the no-slip boundary condition, leads to singular behavior 
of the shearing stress, in the two corners defined by the slider plate and the box walls [1]. Because of the shearing 
motion of the plate, work is done on the fluid, it induces a circulation inside the box, and the fluid is heated as 
a result of viscous dissipation. However, because of the isothermal walls, heat is conducted out of the gas, and a 
steady state is reached after a long time has passed. The velocity field shown in Fig. 10 provides an intuitive physical 
understanding of the flow generated by the sliding plate, which is on top and moves from right to left in the view 
shown. The computational domain was covered by 64x64 square cells, the plate Mach number was set to 1.0, the 
Knudsen number (based on the cavity dimension) was set at 0.005, and the gas was assumed to be air at ambient 
conditions. One of the more interesting results from the solution is the temperature distribution in the box, which 
is shown in Fig. 11. The dimensionless wall temperature is unity in the figure. If no-slip were present, then the 
dimensionless gas temperature would also be unity everywhere along the surface, but it can be seen that a significant 
jump occurs around a great portion of the cavity walls. The same situation develops for the two components of 
velocity. Figure 12 displays the velocity component parallel to the slider plate (near face) which is moving left to right 
in the view shown. If no-slip were present, then the dimensionless gas velocity would be unity (negative) over the 
length of the plate and zero everywhere else on the box boundary. As can be seen, a significant jump in velocity occurs 
over the entire plate and, as a result, it appears to suppress the singular behavior of the stress in the two corners. 
Although the edge-values displayed are actually the values at the midpoints of the cells bordering the walls and one 
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must find the true wall values from extrapolation, the correction is modest and our observations remain unchanged. 
The small ripples seen in the Fig. 12 are believed to be small vortices that are not fully resolved by the grid used. 


VII. CONCLUDING REMARKS 

Our work has been guided by the interest in developing a hybrid method using DSMC and NS, where the focus 
here was limited to a method for handling the NS portion. An Euler scheme introduced by Pullin [15] and further 
developed by Mandal & Deshpande [12] and by Maliett et al. [11], was generalized through use of kinetic-theory 
concepts to cover the NS equations, and second-order solutions were shown to compare well with established second- 
order numerical schemes for the NS equations. Because matching must take place in the near-continuum regime, 
where both DSMC and NS are valid, the NS portion must account for slip at solid boundaries. Use of the kinetic 
split fluxes was shown to lead to a very natural implementation of DSMC type boundary conditions, and these led 
to the appearance of slip. Likewise, use of the kinetic split fluxes at an interface between DSMC and NS in a hybrid 
solution is also expected to be necessary in order to model the correct physics in a nonequilibrium flow, especially 
when determining the inputs to the DSMC portion. In view of this, the equation set (39)-(44), representing KFVS 
for the NS equations, is clearly needed in interfacing with DSMC and in applying boundary conditions at a solid 
surface for the NS portion of a hybrid solution. However, the added step of also using a KFVS scheme for the NS 
solution itself is not an absolute requirement, as one could argue that any numerical solution of the NS equations 
would be acceptable, as long as these kinetic split fluxes were used in applying all boundary conditions. Because the 
DSMC portion is normally expected to take the greater computation time in most problems, one does not need to 
pick a numerical scheme for the NS portion that minimizes time, and therefore, computation time is not an issue. 
This allows one to choose a numerical scheme that offers the greatest compatibility with DSMC, which we believe to 
be the KFVS scheme for the NS equations. On the other hand, in obtaining the data for Figs. 7 and 8, it was found 
that KFVS and Jameson’s SLIP2 required virtually the same computation time. Whether this continues to hold for 
higher spatial dimensions is unknown. 
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Position, x/L 


FIG. 1 . Temperature (T/Tl, - - -) and density ( p/pl , ) ahead of an impulsively started piston in a stationary gas (state 

1 ), where piston Mach number = 5.0, 7 = 5/3 and piston temperature is held fixed. The piston is located at x/L = 1. 



Interface Piston 

FIG. 2. A schematic representation of a hybrid scheme joining the direct simulation Monte Carlo (DSMC) method and a 
Navier-Stokes (NS) numerical method. 
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Mach Number, M 



Mach Number, M 


FIG. 3. Momentum and energy split fluxes for the case of a one-dimensional flow and 7 = 1.4. The momentum flux is 
nondimensionlized by p and the energy flux by p(2RT) * . Comparisons are for the present work, where the assumption 

T^n ~ 0n B ~ — 0.3 was made in Eqs. (40) and (44), ; Mandal-Deshpande KFVS for the Euler equations, - • ; and 

Steger- Warming, . 


L f R 


I 

I 



i— 1 i i + \ i + 1 i + 2 


FIG. 4. Definition of symbols used in applying the MUSCL scheme to the KFVS relations for the NS equations. 
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0.02 



Fluid Velocity, u/ul 


FIG. 5. Normal stress (upper set, r„«/2p) and heat flux (lower set, q n /pc ) in a shock wave versus the fluid velocity, with 
the upstream state designated by 1, A/i = 1.5, and 7 = 1.4. The curves were separated by using a factor of 2 in plotting 

the dimensionless stress. First-order schemes for KFVS, - ■ - Steger- Warming, , and Roe, ••••, are compared with the 

reference solution, . 



FIG. 6. Normal stress (upper set, T nn /2p) and heat flux (lower set, qn/pc) in a shock wave versus the fluid velocity, with 
the upstream state designated by 1, M 1 = 1.5, and 7 = 1.4. The curves were separated by using a factor of 2 in plotting the 
dimensionless stress. Second-order schemes for KFVS, - ■ - •, and Jameson’s SLIP2, - - -, are compared with the reference 
solution, . 


17 







FIG. 7. Shock-tube flow for a pressure ratio of 3, 7 = 1.4, and where the flow is from left to right. Comparisons are for 

second-order KFVS, , Jameson’s SLIP2, , and first-order Roe, . The variables shown are the density p/pl and the 

temperature T/T 1, where state 1 refers to upstream conditions. 




FIG. 8. Shock-tube flow for a pressure ratio of 20, 7 = 1.4, and where the flow is from left to right. Comparisons are for 

second-order KFVS, , Jameson’s SLIP2, , and first-order Roe, ■ • • •. The variables shown are the density p/pl and the 

temperature T/T 1, where state 1 refers to upstream conditions. 
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2 


1.9 

1.8 

1.7 



Position, x/L 

FIG. 9. Temperature ahead of an impulsively started piston in a stationary monatomic gas, where piston Mach number = 
1.0 and where use of the flux boundary condition for an isothermal piston leads to a temperature slip at the piston surface 
(located at x/L — 1). Also shown, by the circle symbol, is Patterson’s theory [11] for temperature slip in a slightly rarefied 
monatomic gas. 



FIG. 10. Velocity field in a square isothermal cavity with a slider plate on top, moving from right to left at a Mach number 
of 1.0. The gas is assumed to be air and Kn (based on the cavity width) = 0.005. 
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X Axis, Cells 


FIG. 11. Temperature distribution in a square isothermal cavity with a slider plate (near face) moving from left to right in 
view shown; conditions same as in Fig. 10. 



X Axis, Cells 


FIG. 12. Velocity component parallel to the slider plate (near face) which is moving from left to right in view shown; 
conditions same as in Fig. 10 with velocity scaled to the sound speed at the wall temperature. 
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APPENDIX B 


A Numerical Study Comparing Kinetic Flux- Vector Splitting for the Navier-Stokes 

Equations with a Particle Method 

T. Lou, D. Dahlby and D. Baganoff 
Department of Aeronautics and Astronautics 
Stanford University, Stanford, CA 94305 
(Submitted to: Journal of Computational Physics, April 30, 1996) 


Numerical solutions based on the method of kinetic flux-vector splitting (KFVS) for the Navier- 
Stokes equations are compared with the direct simulation Monte Carlo method (DSMC) for three 
problems: an impulsively started piston, which emphasizes heat flux; an impulsively started flat 
plate, which emphasizes shearing stress; and a plate sliding past a square cavity, or the lid-driven 
cavity problem, which combines both stress and heat flux. Taking the view that the DSMC method 
provides the correct physical description near material boundaries, the comparisons, which were 
carried out for the conditions of a slightly rarefied flow, show good agreement for temperature and 
velocity slip, and in the prediction of the kinetic split fluxes, verifying the assumptions and the 
approach taken in the KFVS method. 


I. INTRODUCTION 

The theoretical development for the method of kinetic flux-vector splitting (KFVS) for the Navier-Stokes equations 
was introduced in Ref. [4], which represents an extension of work initiated by Deshpande for the Euler equations 
[6,7]. Additionally, it was shown to give good agreement with established numerical schemes for the Navier-Stokes 
equations. Boundary conditions based on the new split fluxes and the kinetic theory were also developed and shown 
to predict slip at a material surface, as a gas becomes rarefied. However, confirmation for the magnitude of the 
predicted slip and the conditions under which the correct predictions are found were not fully explored. Also not 
given was support for the use of a critical approximation, based on the Eucken model, which was introduced to carry 
out the flux-splitting for energy, in the case of a gas having internal structure. The objective of the present work is 
to provide the appropriate analysis by comparing the predictions of the theory presented in [4] with the results of 
simulations carried out with the direct simulation Monte Carlo (DSMC) method [2,3], a method of simulation where 
a large collection of particles is used to model a rarefied gas flow. Of course, the comparisons can only be carried 
out in the near-continuum regime where the computational cost for the DSMC method does not become prohibitive. 
Likewise, for the NS equations to hold, the flow can only be slightly rarefied and the magnitudes of stress and heat 
flux must lie in a range where the occurrence of slip near a solid surface represents the principal modification to the 
fluid physics. However, these conditions are wholly consistent with the objective in [4], where the KFVS method 
was introduced as the continuum counterpart to the DSMC method in an eventual construction of a hybrid scheme 
combining the two. 

The version of the DSMC method used in this study [1,8] divides space into uniform cubical cells. These cells 
are used to identify which particle pairs are candidates for collision during a time step and to compute cell-averaged 
macroscopic quantities at the end of a time step. In general terms, the DSMC method is expected to give reliable 
results when the local mean free path length is large compared with the cell dimension. Particular experience with 
the case of Couette flow has shown that good agreement for viscous stress and heat flux is obtained between DSMC 
and NS when the cell Knudsen number is greater than unity, and progressively more modest agreement is found as 
it is made smaller [5]. In our comparisons, the local cell Knudsen number appearing in the DSMC simulations was 
set at unity or greater for the most dense parts of each flow. The molecular model chosen for the simulations was 
the hard-sphere molecule, for which the transport coefficients vary as the square root of the temperature. In the case 
of a diatomic gas, the vibrational mode was not excited, while the rotational degrees of freedom were set to be in 
equilibrium with the translational degrees of freedom, by setting the so-called collision number to unity. In the DSMC 
simulation,' this leads to the ideal diatomic gas for which 7 = 7/5. 

Two issues are addressed in this study: (i) in the case of a simple gas, the theory in [4] is rather securely founded, 
and therefore, the primary question relates to whether the magnitude of the predicted slip, for the particular flow 
conditions considered, is in agreement with results obtained from DSMC simulations; and (ii) in the case of a gas with 
internal structure, the flux splitting employed in [4] is based on the Eucken approximation, which directly affects the 
predicted split energy fluxes, and the particular approximation used requires confirmation, especially at an isothermal 
material surface where nonequilibrium effects may be quite large. These questions will be investigated by studying 
the highly nonequilibrium flow produced near material surfaces for three problems: an impulsively started piston, 
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which emphasizes heat flux; an impulsively started flat plate, which emphasizes shearing stress; and a plate sliding 
past a square cavity, or the lid-driven cavity problem, which combines both stress and heat flux. 


II. THE KFVS EQUATIONS 

The split kinetic fluxes for the Navier-Stokes equations are given by equations (39)-(44) in Ref. [4] and these are 


reproduced for use here, where the notation employed is the same, and they read: 

F± a „ = pv/rtT/2 [(1 ± a,)S„ ± o 2 (1 - Xi)] (1) 

F*- mm =P (1 ± <*i)(sj + i(l -r„ c „ B )) ± (S„+?£ E ) (2) 

F%_ mm = V2 RT [S,,F± 0 „] + ip [-(1 ± ai )r„ c ,f ± «,«£*] (3) 

Ft-tncm = P'J FT/2 (l±ai) ^S„(i + S 2 ) + X2^ ±Q2(2 + S 2 +X3) (4) 

F? M -. n . rm = + pu„ef„ t ) = i (^) [RTF±„„] (5) 
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S n = u n /V 2ET, S 2 = S 2 „ + 5 t 2 x + Sj 2 


~C£ _ r CE 




fn E = ^/(pVWT) 


= - A' 


CE O) / d\L{ \ 2 (O / 


T ti = p 


dxj dx{ 


A Cartesian coordinates system is assumed in the above equations, with the normal direction represented by n, and the 
two tangential directions by tl and t2. The dimensionless velocity S (speed ratio) determines the two parameters c*i 
and 02) while the dimensionless, Chapman-Enskog expressions for stress t$ e and heat flux q EE are the nonequilibrium 
factors in the quantities xi, X 2 and X3- The sign convention employed assumes the positive split flux F + points in 
the direction of increasing n, and is directed out of a surface enclosing a body of gas. The convention for A - is based 
on the splitting F = F + + F~ , where F is the total flux, and therefore it often evaluates to a negative value. The 
corresponding expression for F^_ mom is not listed as it can be inferred from (3). 

Extreme nonequilibrium conditions in a gas are found near isothermal boundaries, where the one-sided fluxes are 
large, and these provide special conditions for detailed study. The relations developed for an isothermal boundary are 
given by equations (63)-(66) in Ref. [4] and these are also listed below, where again the same notation is followed. 


,) g = Pw V RT w /2tt, 
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(8) 


{Fn-mom) surf ace — Pu//2> 


(Fn — mom)surface — (*S — mom ) 9 5 


(9) 


{F ener gy) sur f ace — {F ener gy)g Pw V ^ FT W / 7 T 


1 + K^r 

4 V 7 ~ 1 


( 10 ) 


In the above relations, the state of the gas near a surface is denoted by subscript g , while the state of the hypothetical 
wall gas is denoted by subscript w. These relations were developed in [4] using DSMC type boundary conditions for 
a nonreacting gas, in which the wall can be viewed as a hypothetical gas at conditions determined by the particular 
boundary conditions employed. In this case, the wall gas is the same gas because of the assumed nonreacting 
interaction between gas molecules. Likewise, the sign convention assumed in the boundary conditions (7)-(10) is one 
where a positive flux points in the direction of positive n, when viewed from the position of the gas at an interface 
with a wall. 


III. IMPULSIVELY STARTED PISTON 


The interest in an impulsively started piston is associated with the fact that the heat transfer rate to an isothermal 
piston can be very high at early time, leading to large nonequilibrium effects. On the other hand, the impulsive 
start requires special attention, as will be seen. Because it is not obvious from the structure of Eqs. (7)-(10) how 
the condition of zero slip is recovered in the continuum limit, and because DSMC simulations become overly costly 
in this limit, it is desirable to further develop the theoretical expressions so that comparisons can be more readily 
assessed. The continuum limit is obtained for conditions of high pressure and large time when the boundary layer is 
relatively thick in relation to the local mean free path length; and for the one-dimensional geometry of an impulsively 
started piston, where the t\ and t2 coordinates are ignored, this leads to small values of the surviving dimensionless 
derivatives t£ e and q % E , with f EE <£ q EE . For calculational purposes, it is easier to consider a moving gas and a 
stationary piston, and therefore, at the piston surface it is appropriate to set the speed ratio to zero, i.e. S n = 0, and 
thus, c*i = 0 and c *2 = I/v^ On using (1) to represent the gas near a surface and on using boundary condition (7), 
the following relations between the conditions in the gas flow and the hypothetical wall-gas values are obtained 



where p = pRT is used for both the gas flow and the hypothetical wall gas. When Eqs. (4)-(6) are substituted into 
boundary condition (10) and on making use of (11) we obtain, for the total energy flux at the surface, the relation 


4(F ener gy) S urf ace _ e~CE < ^ M II j 

Pg^2RTg ~* qn A T g ) 


3 7 - 1 
2\/7r LV 7 - 1 
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The definition for the total energy flux in the coordinate system for which S n = 0 simplifies to the relation 
{Fenergy) surf ace — (fn E , and consequently, the above equation reduces to 





?„ C£ 


37- 1 

7+1 


t CE . 


(13) 


For an isothermal piston, heat conduction represents the dominant effect at large time, and therefore, the above 
equation can be approximated, for f EE <C 1, by 


T g = 1 - (771) t 5 ^"' ' + tf‘) ■ (14) 

which shows that the gas temperature near the surface T g approaches the wall temperature T w , as the magnitudes of the 
nonequilibrium parameters become smaller and smaller in approaching the continuum limit (because of increasing p) 
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Likewise, for positive heat flux (into the wall) the gas temperature is greater than the wall temperature (temperature 
slip). An equivalent slip condition, for the case of slightly rarefied flow of a monatomic gas, was obtained by Patterson 
[9] and it is easy to show that Eq. (14) fully agrees with Patterson’s result (cf. equation (33), p. 125). This agreement 
is found on setting 7 = 5/3, Prandtl number = 2/3, and then approximating (14) for the inverted ratio T g /T w . 
Additionally, it has been shown by Shidlovskiy [10], again for the case of a monatomic gas, that inclusion of the 
thermal accommodation coefficient a introduces a factor (2 — a)/a multiplying the heat flux term (cf. equation (3.16), 
p. 67). Equation (14) was developed here for two purposes: to show that the KFVS formulation (l)-( 6 ) and the 
associated boundary conditions (7)-( 10) agree with related work, and to introduce the extension to the case of a 
polyatomic gas. 

Returning to the more descriptive view, where the gas and piston are both stationary and at one temperature 
before the start of the motion, then the larger the piston velocity after the impulsive start the larger the Mach 
number associated with the shock wave produced, and the greater the changes in density and temperature across 
the gas layer formed near the piston. When one employs the NS (continuum) point of view, then a discontinuity 
appears at the wall in both the temperature and the fluid velocity at time t = 0 + , and the corresponding heat flux 
and normal stress are infinite, and this occurs even for low values of the piston velocity. Clearly, the NS equations 
do not predict the correct physical process at very early time. This raises the interesting question, for this nonsteady 
problem, whether a numerical solution of the NS equations for large time would be independent of developments at 
early time. 

Because Eq. (14) is an analytic result and does not depend explicitly on time, it can be used to qualify a numerical 
solution of the NS equations, because one would expect (14) to provide the correct NS prediction as a numerical 
solution is sequentially improved. The need for a reliable numerical check is the principal reason why the temperature 
ratio in (14) was not inverted, to correspond more directly to Patterson’s expression. On using the variables defined 
by Eq. (14), the analytical relation plots as a straight line, which is shown as a heavy dashed line in Fig. 1 (note: 
q/pc — bq/y/^rj). The four numerical solutions shown in the figure were carried out for the case of a monatomic gas 
and a piston Mach number of unity ( M s h oc k — 1 869), using a second-order finite-volume scheme (first-order time) 
together with a range of cell sizes (see Ref. [4] for identification of the scheme used). Time appears as a parameter 
along each curve, with large time corresponding to small values of the abscissa. The physical scale is set by the values 
of the undisturbed density po, speed of sound Co, coefficient of viscosity po, and cell length Ax. These quantities can be 
used to either define a reference cell Reynolds number or a reference cell Knudsen number through the kinetic-theory, 
hard-sphere relation 

(15) 


where C = \/SRT / 7 r is the mean thermal speed and A is the mean free path length. Because we are interested in 
comparisons with DSMC, it is physically more meaningful to use a reference cell Knudsen number, Kuq = Ao/Ax, 
with x = x/\q , t = tCo/Xo , and 



At — AtCo/ A xKtiq — 


—{CFL/I<n) o, 

7T7 


(16) 


as measures of the physical scale. In Fig. 1 the numerical solutions are shown for Ktiq = 2, 4, 8, 16 and we see that the 
largest value is needed to get good agreement with (14). Because the gas density near the piston surface is nearly four 
times the undisturbed density (see Fig. 3), this translates into Kn wa u « 4 and therefore we must have Ax < X wa ii /4 
to obtain a reliable NS solution at early time. The conditions near t = 0 + necessitated the use of a very small Af. 
This led to the use of values of CFL ranging from 0.15 to 0.0375 as the Knudsen number was increased. At this point 
we do not know whether the range of the independent variable displayed in Fig. 1 corresponds to conditions where the 
NS system predicts the correct physics, only that the values employed are required for a consistent numerical solution 
of the NS equations. 

In turning to a DSMC simulation, it is clear that the same reference cell Knudsen number, /in 0 = 16, should 
initially be used in making a comparison. Figure 2 shows such a comparison for two times: a time at which the shock 
wave and the thermal layer are both still forming, t = 4.64 (1,600 time steps NS, 400 DSMC); and the time at which 
they just begin to separate, t — 11.6 (4,000 time steps NS, 1,000 DSMC). Numerical instability with NS at early time 
necessitated a smaller time step (factor of 4) than that used with DSMC. The dimensionless time employed is based 
on the collision time in the reference state and so these times are truly short. At the early time the temperature 
profiles match somewhat poorly overall, while at the later time the thermal layers, alone, begin to match rather well. 
Because density is a less sensitive variable, the match at both times is surprisingly good, considering the extremely 
short time represented. Because the DSMC method is computationally intensive, it was necessary to increase the 
values of Af and Ax by a factor of 4, and decrease the value of I\ no by the same factor, to study still larger time. 
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Figure 3 makes the same comparison after the shock wave and the thermal layer have clearly separated, t = 46.4 
(16,000 time steps NS; 1,000 DSMC), and it is seen that agreement is very good, except at the shock front itself, for 
which it is well known that NS gives a poor prediction for shock- wave profile. 

The inherent statistical fluctuations which are characteristic of the DSMC method, especially for a nonsteady 
problem for which extended time averaging is not possible, do not allow for a detailed study of small differences 
represented by the temperature slip seen in Figs. 2 and 3. Nevertheless, in Fig. 4 we attempt to make a comparison of 
the time dependent temperature slip at the piston surface for the two methods. In this case, the DSMC results were 
time averaged over a small local interval about each plotted point to reduce statistical scatter. These simulations 
were carried out with a number density of approximately 8,000 particle per cell near the piston surface and roughly 
one million particles for the entire simulation. DSMC results were obtained for Kn o = 4,8,16,32 by sequentially 
reducing the cell size by a factor of 2 while holding A 0 fixed. These data include the runs shown in Figs. 2 and 3, and 
likewise, are for the same conditions. Because of the varying cell size, it would be necessary to extrapolate the data 
for each run to the position of the piston surface in order to produce a consistent display, but this approach amplifies 
statistical scatter and proves impractical for such a time-dependent simulation with DSMC. The alternative was to 
select the center position of the largest DSMC cell (Ktiq — 4) as the reference and display the data for all runs for 
that position. This approach has the added feature that it provides a consistency check on the DSMC method itself. 
As seen in Fig. 4, the DSMC data overlap nicely, demonstrating that convergence has been obtained. Two curves are 
shown for the NS calculation: the dashed curve was obtained by extrapolating the data to the piston surface; and 
the solid curve represents the value at the center position of the largest cell used in the DSMC method (A no = 4). 
For the NS calculation, the cell size for I<n 0 = 16 was used. The separation between the two curves shows that the 
temperature gradient near the piston surface is very steep necessitating the use of a small cell size. In view of the 
different curves displayed, the solid curve for NS should be compared with the DSMC data, which appears to suggests 
that the boundary conditions (7)-(ll) for NS over predict the temperature slip for these conditions. For large time 
the two should agree fully, but it did not appear feasible to extend the DSMC runs to verify this with the computer 
workstation employed in this part of the study. When considering the fact that the NS equations are not expected 
to represent the correct physics for large nonequilibrum (for example, q/pc > 0.1 for t < 50, see Figs. 1 and 4), the 
agreement seen in Figs. 2-4 is very encouraging, since it confirms the accepted view that slip conditions at a surface 
represent the principal corrections needed to be added to the NS system when dealing with a slightly rarefied flow. 

If any discrepancy exists between the values of the split kinetic fluxes defined by Eqs. (l)-(6) and the corresponding 
values from DSMC simulations, then the differences should be seen at the piston surface where nonequilibrium is 
the greatest. Fig. 5 presents the corresponding comparisons for the mass and energy split fluxes, and shows that 
the agreement, for a monatomic gas, is extremely good. In the DSMC simulations the positive split fluxes were 
obtained by monitoring the passage of individual particles as they left the gas and crossed the piston surface, while 
the negative split fluxes were obtained by monitoring the particle emission from the piston surface introduced by 
the DSMC boundary conditions. For the NS solution, the state of the flow from the numerical solution was used 
to evaluate the positive split fluxes using the defining equations (l)-(6). Similarly, the conditions representing the 
isothermal piston were used to compute the negative split fluxes. 

In the case of a polyatomic gas, the Eucken approximation was introduced in [4] to develop the split fluxes for 
energy; and it is of interest to determine whether the particular approximation used is supported by DSMC. The 
same simulations were repeated for the case of an ideal diatomic gas, assuming rotational degrees of freedom are in 
equilibrium with translation (7 = 7/5), and a comparison of split fluxes for momentum and energy is shown in Fig. 6. 
The excellent agreement seen confirms that the Eucken approximation, as implemented in [4], is capturing the proper 
physics in the splitting of energy flux, for both the translational and rotational components. 


IV. IMPULSIVELY STARTED FLAT PLATE 

Viscous stress becomes the dominant nonequilibrium effect for an impulsively started flat plate, which provides 
an alternate environment for comparison. For the one-dimensional geometry of an infinite, impulsively started flat 
isothermal plate moving parallel to its surface in, say, the t\ direction, the t\ and £2 coordinates may be disregarded 
in (l)-(6); and the shearing stress f„ t f , along with the normal heat-flux component q^ E , become the principal 
nonequilibrium quantities. Again, for calculational purposes it is easier to consider a moving gas and a stationary 
plate. At the plate surface we may then set S n = 0 and thus a\ = 0 and — l/y/n- On using Eq. (1) to represent 
the gas near the plate and on using boundary condition (7), exactly the same density and temperature relations are 
found as for the impulsively started piston, given by (11). 

On using (1) and (3) to represent the gas near a stationary surface, in which we set qf x E to zero because of uniformity, 
and on using boundary condition (9), we find 
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The definition for the total flux of transverse momentum leads to the substitution (F t \-mom) surf ace — — TVui, because 
S n = 0. On solving for S t i, we then obtain the relation 


Sti 



(18) 


The quantity S t i represents the dimensionless velocity of the fluid next to a stationary surface, i.e., velocity slip. 
Although Eq. (18) does not depend directly on the value of the ratio of specific heats 7 , if the velocity slip is based on 
the reference plate speed u p i atey and a plate Mach number M p ute is introduced, then the factor y/rr/2 is replaced by 
y/n/2~f /M p i ate and the 7 dependence becomes evident. Here again, if we neglect r£ E then the above relation can be 
shown to agree fully with Patterson’s result [9] for velocity slip (cf. equation (31), p. 125). If boundary condition (10) 
is handled in the same way as for the piston, and if the quadratic terms Stxf nt i and 5 t 2 x are dropped, then exactly 
the same relation for the temperature slip as for the impulsively started piston is gotten, i.e., Eq. (14). 

Just as for the case of the impulsively started piston, Eq. (18) can be used to qualify a numerical solution of the 
NS equations for an impulsively started flat plate. Except for a greater sensitivity to the impulsive start, which 
requires the use of still smaller time steps as the reference cell Knudsen number is increased, the conclusions drawn 
are essentially the same as those found in the study that led to the data presented in Fig. 1, and will not be repeated. 
Likewise, in order to emphasize nonequilibrium effects, an isothermal plate and a Mach number of unity were selected 
as appropriate conditions for study. 

Using the more descriptive frame of reference where the gas is initially stationary and the plate is given an impulsive 
start, NS and DSMC results for velocity are displayed and compared in Fig. 7, at the dimensionless times t = 12.4 
and t — 68.0. As can be seen, even though the velocity slip is fairly large at these short times, the NS solution 
compares extremely well with the DSMC results. In Fig. 8 the gas velocity at the surface of the plate is displayed as 
a function of time, showing that it approaches the plate velocity asymptotically. The DSMC results were obtained 
for Kno = 4 and 8 ; and data for both runs were plotted for the location corresponding to the center position of the 
largest cell, Kno = 4. The complete overlap of the symbols clearly demonstrates consistency, or convergence, in the 
DSMC results. In the case of the NS solution, the data for I<n 0 = 12 were used and extrapolated to the plate surface 
(dashed curve); the data were also evaluated at the center position of the largest cell (solid curve) used in the DSMC 
simulations. As the gradient near the plate is more modest here, the agreement is sufficiently close that one does 
not have to distinguish between the different curves. Equation (18) shows that the magnitude of f£ E is virtually the 
same as the slip (1 - u/u p i ate ) seen in the figure. Therefore, the degree of nonequilibirum is quite large, yet the NS 
prediction for velocity slip agrees very nicely with DSMC for these rather extreme conditions. 

Continuing with comparisons for the split kinetic fluxes evaluated at the plate surface, Fig. 9 presents the results for 
the case of a monatomic gas, and the agreement seen is remarkably good. Likewise, Fig. 10 presents the split fluxes 
for a diatomic gas, showing equally good agreement. These quantities are displayed for a frame of reference where 
the gas is initially stationary and the plate is given an impulsive start. Therefore, asymptotic results deduced from 
Eq. (l)-( 6 ) must be transformed to obtain the limiting values seen in Figs. 9 and 10. Beyond the excellent comparisons 
seen, the most important observation relates to the component quantities making up the energy split fluxes for the 
diatomic gas. It is clear that the translational and rotational degrees of freedom are being properly handled, and 
therefore, the approach used in Ref. [4], in introducing the Eucken approximation, appears to be working well. 


V. LID-DRIVEN CAVITY FLOW 

In the two cases studied above, heat flux and viscous stress were separately dominant; but both can become 
important in the lid-driven cavity problem, and an element of complexity is added by the two dimensions of the flow. 
However, focus will be placed on the steady-state condition for which the CFL restriction should not be as severe as 
for the case of an impulsive start. Here again, it is useful to assume isothermal surfaces to produce a large degree 
of nonequilibrium. On the other hand, the lid Mach number was set to Mud = 0.5 so that nonequilibrium effects 
in the corners for steady state would not be unduly large. The reference Knudsen number, based on the dimension 
of the square cavity L, and defined by Kn L = \q/L, was set equal to 0.01, where Ao is the mean free path length 
evaluated at the wall temperature and the initial state of the gas. The value of the Knudsen number was chosen so 
as to correspond to the near-continuum regime where NS is expected to be valid and the DSMC simulation does not 
become overly intensive. In the following, discussion will be limited to the case of a diatomic gas, as the monatomic 
case has been adequately covered above. 
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A steady state, two-dimensional NS solution, based on a 128 x 128 square mesh, is presented in Fig. 11 showing the 
component of velocity lying parallel to the lid. In the figure, the lid is on the near face and moves from left to right, 
for which the velocity is defined to be negative. The effect of velocity slip is clearly seen, both on the lid itself and at 
the two corners formed by the lid and walls. No slip would correspond to the magnitude of the dimensionless fluid 
velocity u/cq being equal to the lid Mach number, in this case 0.5. The view shown is useful in serving as a mental 
aid in presenting the comparisons to be reviewed below. For example, in the views that follow, the lid together with 
the two faces on either side will be unwrapped and displayed in planar form when various boundary quantities are 
compared. 

In order to judge the validity of the numerical solution presented in Fig. 11, analytical relations similar to Eqs. (14) 
and (18) are needed. Using a coordinate system where n is taken to be perpendicular to the lid, tl is parallel to 
the lid and t2 ignorable, then at various points along the lid one would expect the stresses t ee and r EE and the 
heat flux components q EE and to be important. On this basis the slip relations found above may have to be 
generalized. When the algebra leading to (11) is repeated using S n — 0 alone, exactly the same relations are found 
A slight generalization to (18) is required, given by 
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which introduces the heat flux component aligned with the lid, a quantity that may be important in the corners. 
More terms must be retained in the generalization of (13), for the temperature slip, which becomes 
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In developing Eq. (20), the coordinate system in which the gas is moving and the wall is stationary was again used. 
Therefore, a transformation must be introduced when (20) is used to analyze the moving lid. Because of the complexity 
of (19) and (20), the approach employed in Fig. 1 is less useful here. Here too, Eqs. (19) and (20) can be shown to 
reduce to Patterson’s results [9]. 

In the coordinate system defined by Fig. 11, the velocity component u lies parallel to the lid while the velocity 
component v lies parallel to the two faces on either side of the lid. Therefore, if we unwrap the adjacent faces and 
display the tangential velocity along the surface as a function of the surface position for the NS solution, we obtain 
the function seen in the top view of Fig. 12, shown as a solid curve. This function contains both u and v and is not 
merely a copy of the edge values for u shown in Fig. 11. The corresponding theoretical prediction is given by Eq. (19) 
and is shown as the heavy dotted curve (a transformation must be applied to (19) to obtain the values along the lid). 
It is clear that agreement is only found outside the two corner regions. This can be understood by reviewing the plot 
shown in the bottom view of the figure, which gives the normal velocity component 5 n as a function of s. Equation 
(19) was derived on the basis of the assumption S n = 0 and it is clear from the plot for S n that the assumption does 
not hold in the two corners. Therefore, use of Eq. (19) as a check on the validity of the numerical solution obtained 
must be confined to areas near s = 0.5, 1.5, 2.5; and this check clearly shows that reliable numerical results were in 
fact obtained. 

Practical considerations made it necessary to limit the DSMC simulation to a 64 x 64 array of cells. Past experience 
with the DSMC method for a steady state problem led to the decision to use an average number density of approx- 
imately 64 particles per cell, leading to a total of roughly 0.25 x 10 6 particles employed in the simulation. Roughly 
12,000 time steps were used in the time averaging of the data which gave a sample size for each cell of approximately 
0.75 x 10 6 . In this case it was more appropriate to run the simulation on a parallel supercomputer, which required a 
total of 116 node-hours to carry out the comparisons. The top view of Fig. 13 compares the NS solution against DSMC 
results for the tangential velocity versus the surface position. In each case the data were projected to the position of 
the surface for comparison. As can be seen, the correspondence is quite complete, even including sharp spikes in the 
two corners. Consequently, the NS prediction for the slip velocity is quite outstanding for these conditions. Because of 
the low Mach number chosen, Mud = 0.5, the temperature rise is fairly small (less than 5%) and the DSMC data for 
temperature exhibit considerable statistical scatter, as is seen in the bottom view in the figure, and a fully equivalent 
judgment concerning the NS temperature slip, DSMC results, and Eq. (20) cannot be made. 

Because the wall temperature is specified for an isothermal wall, the NS solution only controls the density of the 
hypothetical wall gas p w (see Eq. (11)). Furthermore, because the emission from the wall is controlled by a Maxwellian 
distribution in the particular application (isothermal wall) of the KFVS method being considered, the split fluxes 
directed out of the wall are not overly sensitive to the NS solution. Thus, it makes sense to focus attention on the 
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split fluxes directed out of the gas and into the wall. Figure 14 gives a comparison for this single component of the 
two kinetic split fluxes for energy, where the upper set is for the translational energy component and the lower set is 
for the rotational energy component. As seen, the NS predictions compare extremely well with the DSMC results, 
including the excursions in the two corners. Likewise, the good match for rotational energy shows that the form of the 
Eucken approximation employed in the development of the KFVS method proves to be valid even for these conditions. 

Although prediction of the values of the fluid variables at the lid and walls represents a more severe test of the 
KFVS theory than that for the interior regions of the lid-driven cavity problem, it is also of interest to consider 
one comparison for the entire flow. Figure 15 displays the v-component of velocity for a number of transverse slices 
positioned along the x axis. The NS solution is given by the solid curves and the DSMC results by the symbols, 
likewise showing good agreement between the two. 


VI. CONCLUDING REMARKS 

The method of kinetic flux-vector splitting for the Navier-Stokes equations was introduced primarily as the contin- 
uum counterpart to the DSMC method in the eventual development of a hybrid scheme [4], A principal requirement 
in joining the two methods at a fluid interface is the presence of compatible split fluxes. Because NS is not valid 
in rarefied flow and the use of the DSMC method becomes overly costly in the deep continuum, matching must be 
carried out in the near-continuum where the flow is only slightly rarefied, a degree of rarefaction which may be defined 
as the regime where the local cell Knudsen number is of order unity. These conditions were considered in selecting 
the test conditions reported above and very good agreement was found for the split fluxes in the two schemes, even 
for the extreme nonequilibrium conditions found near isothermal surfaces, conditions surely more severe than those 
found at most any other interface located within the flow. Because the split fluxes for mass and momentum do not 
depend on 7 , it is certainly expected that one would obtain equally good results for monatomic and diatomic gasses 
for these quantities. However, the split fluxes for energy clearly depend on the additional internal energy carried by 
polyatomic molecules, and for this case it was necessary to make use of a particular interpretation of the Eucken 
model to carry out the splitting (see Ref. [4]). Therefore, comparisons such as those seen in Figs. 10 and 14 prove 
to be of great value in justifying the assumptions made. Beyond the fact that the split fluxes defined in KFVS and 
DSMC have been shown to agree remarkable well, which is an important step in the development of a hybrid scheme, 
it was also shown that the slip conditions for temperature and velocity at a material surface also agreed quite well. 
This is also a significant result because the flow is fully expected to be slightly rarefied near a material surface for 
most conditions for which a hybrid scheme would be employed. 
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FIG. 1. Gas temperature at the surface of an impulsively started isothermal piston moving into a stationary monatomic 

gas. Theoretical relation given by Eq. (14), . Numerical integration of the NS equations for M v i 3t0 n — 1, 7 = 5/3, and 

reference cell Knudsen number Kno — 2, . . . . ; 4, ; 8, - . - . - ; and 16, . 




FIG. 2. Gas temperature and density ahead of an isothermal piston (located at x = 0) for NS (solid curves) and DSMC 
(symbols), at two early times in the formation of the shock wave and the thermal layer, for M ptst on = 1 and 7 = 5/3. 
Dimensionless times correspond to t = 4.64 and t = 11.6. 
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FIG. 3. Gas temperature and density ahead of an isothermal piston (located at x = 0) for NS (solid curves) and DSMC 
(symbols), for M pia ton = 1 , 7 = 5/3, and at a time t = 46.4 when the shock wave and thermal layer have clearly separated, 
showing a well defined thermal layer. The DSMC results are for Kno = 4 while the NS results are for Kno — 16, i.e., Ax for 
NS is 4 times smaller than for DSMC. 



FIG. 4. Gas temperature at the surface of an isothermal piston versus dimensionless time, for M p i st on — 1 and 7 = 5/3. 
DSMC cell dimensions correspond to Kno = 2, o; 4, x ; 8, *; 16, +. The NS solution was evaluated for Kno = 16 and 
projected to the piston surface (dashed curve); same solution was evaluated at cell center for the largest DSMC cell (solid 
curve). 
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FIG. 5. Kinetic split fluxes evaluated at the surface of an isothermal piston, for Mpiaton — 1 and 7 = 5/3. The numerical 
solution of the NS equations (solid curves) tire compared with results from a DSMC simulation (symbols). Mass and energy 
split fluxes are referenced to the values poCo and PoCq , respectively. 
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FIG. 6. Kinetic split fluxes evaluated at the surface of an isothermal piston, for M pt3 ton = 1 and 7 = 7/5. Momentum and 
energy split fluxes are referenced to the values pocl and pqCq, respectively. 
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FIG. 7. Velocity distribution above an impulsively started isothermal flat plate, for M p i a te — 1 and 7 = 5/3. The two times 
shown cire for t = 12.4 and 68.0, with NS solution (solid curves) and DSMC simulation (symbols). 



FIG. 8. Gas velocity at the surface of an impulsively started isothermal flat plate versus time, for M p i a te = 1 and 7 = 5/3. 
NS solution for Kno — 12, with extrapolation to surface (dashed curve) and location at cell center for largest DSMC cell (solid 
curve). DSMC simulation for Kno =4, o ; 8, +• 
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FIG. 9. Kinetic split fluxes evaluated at the surface of an isothermal flat plate, for M p i a te = 1 and 7 = 5/3. The numerical 
solution of the NS equations (solid curves) are compared with results from a DSMC simulation (symbols). Momentum and 
energy split fluxes are referenced to the values PqCq and po Cq , respectively. 
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FIG. 10. Kinetic split fluxes evaluated at the surface of an isothermal flat plate, for M p /ote = 1 and 7 = 7/5. The numerical 
solution of the NS equations (solid curves) are compared with results from a DSMC simulation (symbols). Momentum and 
energy split fluxes are referenced to the values poc% and poCo, respectively. 
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FIG. 11 . Variation of the u-component of velocity for the NS solution to the lid-driven cavity problem. Lid velocity is from 
left to right on the near face, with Mud = 0.5 and 7 = 7/5. 




FIG. 12. Tangential and normal velocities at the surface of the lid and cavity for the NS solution. Upper view: numerical 
solution of the NS equations (solid curve) and the theoretical expression, Eq. (19), for the tangential velocity tq/co (dotted 
curve). Lower view: velocity component normal to the surface of the cavity, S„, for the NS solution. 
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FIG. 13. Tangential velocity, u t /co t and temperature, T/To , at the surface of the lid and cavity for the NS solution (solid 
curves) versus result for a DSMC simulation (symbols). 



FIG. 14. Translational (upper curves) and rotational (lower curves) energy split fluxes at the surface of the lid and cavity 
for NS (solid curves) and DSMC (symbols). Energy split fluxes are referenced to the value PqCq. 
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FIG. 15. Comparison of the v-component of velocity, v/cq , in the cavity for NS (solid curves) and DSMC (symbols). 
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